1 use num::{One, Zero}; 2 #[cfg(feature = "abomonation-serialize")] 3 use std::io::{Result as IOResult, Write}; 4 5 use approx::{AbsDiffEq, RelativeEq, UlpsEq}; 6 use std::any::TypeId; 7 use std::cmp::Ordering; 8 use std::fmt; 9 use std::hash::{Hash, Hasher}; 10 use std::marker::PhantomData; 11 use std::mem; 12 13 #[cfg(feature = "serde-serialize-no-std")] 14 use serde::{Deserialize, Deserializer, Serialize, Serializer}; 15 16 #[cfg(feature = "abomonation-serialize")] 17 use abomonation::Abomonation; 18 19 use simba::scalar::{ClosedAdd, ClosedMul, ClosedSub, Field, SupersetOf}; 20 use simba::simd::SimdPartialOrd; 21 22 use crate::base::allocator::{Allocator, SameShapeAllocator, SameShapeC, SameShapeR}; 23 use crate::base::constraint::{DimEq, SameNumberOfColumns, SameNumberOfRows, ShapeConstraint}; 24 use crate::base::dimension::{Dim, DimAdd, DimSum, IsNotStaticOne, U1, U2, U3}; 25 use crate::base::iter::{ 26 ColumnIter, ColumnIterMut, MatrixIter, MatrixIterMut, RowIter, RowIterMut, 27 }; 28 use crate::base::storage::{Owned, RawStorage, RawStorageMut, SameShapeStorage}; 29 use crate::base::{Const, DefaultAllocator, OMatrix, OVector, Scalar, Unit}; 30 use crate::{ArrayStorage, SMatrix, SimdComplexField, Storage, UninitMatrix}; 31 32 use crate::storage::IsContiguous; 33 use crate::uninit::{Init, InitStatus, Uninit}; 34 #[cfg(any(feature = "std", feature = "alloc"))] 35 use crate::{DMatrix, DVector, Dynamic, VecStorage}; 36 use std::mem::MaybeUninit; 37 38 /// A square matrix. 39 pub type SquareMatrix<T, D, S> = Matrix<T, D, D, S>; 40 41 /// A matrix with one column and `D` rows. 42 pub type Vector<T, D, S> = Matrix<T, D, U1, S>; 43 44 /// A matrix with one row and `D` columns . 45 pub type RowVector<T, D, S> = Matrix<T, U1, D, S>; 46 47 /// The type of the result of a matrix sum. 48 pub type MatrixSum<T, R1, C1, R2, C2> = 49 Matrix<T, SameShapeR<R1, R2>, SameShapeC<C1, C2>, SameShapeStorage<T, R1, C1, R2, C2>>; 50 51 /// The type of the result of a matrix sum. 52 pub type VectorSum<T, R1, R2> = 53 Matrix<T, SameShapeR<R1, R2>, U1, SameShapeStorage<T, R1, U1, R2, U1>>; 54 55 /// The type of the result of a matrix cross product. 56 pub type MatrixCross<T, R1, C1, R2, C2> = 57 Matrix<T, SameShapeR<R1, R2>, SameShapeC<C1, C2>, SameShapeStorage<T, R1, C1, R2, C2>>; 58 59 /// The most generic column-major matrix (and vector) type. 60 /// 61 /// # Methods summary 62 /// Because `Matrix` is the most generic types used as a common representation of all matrices and 63 /// vectors of **nalgebra** this documentation page contains every single matrix/vector-related 64 /// method. In order to make browsing this page simpler, the next subsections contain direct links 65 /// to groups of methods related to a specific topic. 66 /// 67 /// #### Vector and matrix construction 68 /// - [Constructors of statically-sized vectors or statically-sized matrices](#constructors-of-statically-sized-vectors-or-statically-sized-matrices) 69 /// (`Vector3`, `Matrix3x6`…) 70 /// - [Constructors of fully dynamic matrices](#constructors-of-fully-dynamic-matrices) (`DMatrix`) 71 /// - [Constructors of dynamic vectors and matrices with a dynamic number of rows](#constructors-of-dynamic-vectors-and-matrices-with-a-dynamic-number-of-rows) 72 /// (`DVector`, `MatrixXx3`…) 73 /// - [Constructors of matrices with a dynamic number of columns](#constructors-of-matrices-with-a-dynamic-number-of-columns) 74 /// (`Matrix2xX`…) 75 /// - [Generic constructors](#generic-constructors) 76 /// (For code generic wrt. the vectors or matrices dimensions.) 77 /// 78 /// #### Computer graphics utilities for transformations 79 /// - [2D transformations as a Matrix3 <span style="float:right;">`new_rotation`…</span>](#2d-transformations-as-a-matrix3) 80 /// - [3D transformations as a Matrix4 <span style="float:right;">`new_rotation`, `new_perspective`, `look_at_rh`…</span>](#3d-transformations-as-a-matrix4) 81 /// - [Translation and scaling in any dimension <span style="float:right;">`new_scaling`, `new_translation`…</span>](#translation-and-scaling-in-any-dimension) 82 /// - [Append/prepend translation and scaling <span style="float:right;">`append_scaling`, `prepend_translation_mut`…</span>](#appendprepend-translation-and-scaling) 83 /// - [Transformation of vectors and points <span style="float:right;">`transform_vector`, `transform_point`…</span>](#transformation-of-vectors-and-points) 84 /// 85 /// #### Common math operations 86 /// - [Componentwise operations <span style="float:right;">`component_mul`, `component_div`, `inf`…</span>](#componentwise-operations) 87 /// - [Special multiplications <span style="float:right;">`tr_mul`, `ad_mul`, `kronecker`…</span>](#special-multiplications) 88 /// - [Dot/scalar product <span style="float:right;">`dot`, `dotc`, `tr_dot`…</span>](#dotscalar-product) 89 /// - [Cross product <span style="float:right;">`cross`, `perp`…</span>](#cross-product) 90 /// - [Magnitude and norms <span style="float:right;">`norm`, `normalize`, `metric_distance`…</span>](#magnitude-and-norms) 91 /// - [In-place normalization <span style="float:right;">`normalize_mut`, `try_normalize_mut`…</span>](#in-place-normalization) 92 /// - [Interpolation <span style="float:right;">`lerp`, `slerp`…</span>](#interpolation) 93 /// - [BLAS functions <span style="float:right;">`gemv`, `gemm`, `syger`…</span>](#blas-functions) 94 /// - [Swizzling <span style="float:right;">`xx`, `yxz`…</span>](#swizzling) 95 /// 96 /// #### Statistics 97 /// - [Common operations <span style="float:right;">`row_sum`, `column_mean`, `variance`…</span>](#common-statistics-operations) 98 /// - [Find the min and max components <span style="float:right;">`min`, `max`, `amin`, `amax`, `camin`, `cmax`…</span>](#find-the-min-and-max-components) 99 /// - [Find the min and max components (vector-specific methods) <span style="float:right;">`argmin`, `argmax`, `icamin`, `icamax`…</span>](#find-the-min-and-max-components-vector-specific-methods) 100 /// 101 /// #### Iteration, map, and fold 102 /// - [Iteration on components, rows, and columns <span style="float:right;">`iter`, `column_iter`…</span>](#iteration-on-components-rows-and-columns) 103 /// - [Elementwise mapping and folding <span style="float:right;">`map`, `fold`, `zip_map`…</span>](#elementwise-mapping-and-folding) 104 /// - [Folding or columns and rows <span style="float:right;">`compress_rows`, `compress_columns`…</span>](#folding-on-columns-and-rows) 105 /// 106 /// #### Vector and matrix slicing 107 /// - [Creating matrix slices from `&[T]` <span style="float:right;">`from_slice`, `from_slice_with_strides`…</span>](#creating-matrix-slices-from-t) 108 /// - [Creating mutable matrix slices from `&mut [T]` <span style="float:right;">`from_slice_mut`, `from_slice_with_strides_mut`…</span>](#creating-mutable-matrix-slices-from-mut-t) 109 /// - [Slicing based on index and length <span style="float:right;">`row`, `columns`, `slice`…</span>](#slicing-based-on-index-and-length) 110 /// - [Mutable slicing based on index and length <span style="float:right;">`row_mut`, `columns_mut`, `slice_mut`…</span>](#mutable-slicing-based-on-index-and-length) 111 /// - [Slicing based on ranges <span style="float:right;">`rows_range`, `columns_range`…</span>](#slicing-based-on-ranges) 112 /// - [Mutable slicing based on ranges <span style="float:right;">`rows_range_mut`, `columns_range_mut`…</span>](#mutable-slicing-based-on-ranges) 113 /// 114 /// #### In-place modification of a single matrix or vector 115 /// - [In-place filling <span style="float:right;">`fill`, `fill_diagonal`, `fill_with_identity`…</span>](#in-place-filling) 116 /// - [In-place swapping <span style="float:right;">`swap`, `swap_columns`…</span>](#in-place-swapping) 117 /// - [Set rows, columns, and diagonal <span style="float:right;">`set_column`, `set_diagonal`…</span>](#set-rows-columns-and-diagonal) 118 /// 119 /// #### Vector and matrix size modification 120 /// - [Rows and columns insertion <span style="float:right;">`insert_row`, `insert_column`…</span>](#rows-and-columns-insertion) 121 /// - [Rows and columns removal <span style="float:right;">`remove_row`, `remove column`…</span>](#rows-and-columns-removal) 122 /// - [Rows and columns extraction <span style="float:right;">`select_rows`, `select_columns`…</span>](#rows-and-columns-extraction) 123 /// - [Resizing and reshaping <span style="float:right;">`resize`, `reshape_generic`…</span>](#resizing-and-reshaping) 124 /// - [In-place resizing <span style="float:right;">`resize_mut`, `resize_vertically_mut`…</span>](#in-place-resizing) 125 /// 126 /// #### Matrix decomposition 127 /// - [Rectangular matrix decomposition <span style="float:right;">`qr`, `lu`, `svd`…</span>](#rectangular-matrix-decomposition) 128 /// - [Square matrix decomposition <span style="float:right;">`cholesky`, `symmetric_eigen`…</span>](#square-matrix-decomposition) 129 /// 130 /// #### Vector basis computation 131 /// - [Basis and orthogonalization <span style="float:right;">`orthonormal_subspace_basis`, `orthonormalize`…</span>](#basis-and-orthogonalization) 132 /// 133 /// # Type parameters 134 /// The generic `Matrix` type has four type parameters: 135 /// - `T`: for the matrix components scalar type. 136 /// - `R`: for the matrix number of rows. 137 /// - `C`: for the matrix number of columns. 138 /// - `S`: for the matrix data storage, i.e., the buffer that actually contains the matrix 139 /// components. 140 /// 141 /// The matrix dimensions parameters `R` and `C` can either be: 142 /// - type-level unsigned integer constants (e.g. `U1`, `U124`) from the `nalgebra::` root module. 143 /// All numbers from 0 to 127 are defined that way. 144 /// - type-level unsigned integer constants (e.g. `U1024`, `U10000`) from the `typenum::` crate. 145 /// Using those, you will not get error messages as nice as for numbers smaller than 128 defined on 146 /// the `nalgebra::` module. 147 /// - the special value `Dynamic` from the `nalgebra::` root module. This indicates that the 148 /// specified dimension is not known at compile-time. Note that this will generally imply that the 149 /// matrix data storage `S` performs a dynamic allocation and contains extra metadata for the 150 /// matrix shape. 151 /// 152 /// Note that mixing `Dynamic` with type-level unsigned integers is allowed. Actually, a 153 /// dynamically-sized column vector should be represented as a `Matrix<T, Dynamic, U1, S>` (given 154 /// some concrete types for `T` and a compatible data storage type `S`). 155 #[repr(C)] 156 #[derive(Clone, Copy)] 157 pub struct Matrix<T, R, C, S> { 158 /// The data storage that contains all the matrix components. Disappointed? 159 /// 160 /// Well, if you came here to see how you can access the matrix components, 161 /// you may be in luck: you can access the individual components of all vectors with compile-time 162 /// dimensions <= 6 using field notation like this: 163 /// `vec.x`, `vec.y`, `vec.z`, `vec.w`, `vec.a`, `vec.b`. Reference and assignation work too: 164 /// ``` 165 /// # use nalgebra::Vector3; 166 /// let mut vec = Vector3::new(1.0, 2.0, 3.0); 167 /// vec.x = 10.0; 168 /// vec.y += 30.0; 169 /// assert_eq!(vec.x, 10.0); 170 /// assert_eq!(vec.y + 100.0, 132.0); 171 /// ``` 172 /// Similarly, for matrices with compile-time dimensions <= 6, you can use field notation 173 /// like this: `mat.m11`, `mat.m42`, etc. The first digit identifies the row to address 174 /// and the second digit identifies the column to address. So `mat.m13` identifies the component 175 /// at the first row and third column (note that the count of rows and columns start at 1 instead 176 /// of 0 here. This is so we match the mathematical notation). 177 /// 178 /// For all matrices and vectors, independently from their size, individual components can 179 /// be accessed and modified using indexing: `vec[20]`, `mat[(20, 19)]`. Here the indexing 180 /// starts at 0 as you would expect. 181 pub data: S, 182 183 // NOTE: the fact that this field is private is important because 184 // this prevents the user from constructing a matrix with 185 // dimensions R, C that don't match the dimension of the 186 // storage S. Instead they have to use the unsafe function 187 // from_data_statically_unchecked. 188 // Note that it would probably make sense to just have 189 // the type `Matrix<S>`, and have `T, R, C` be associated-types 190 // of the `RawStorage` trait. However, because we don't have 191 // specialization, this is not bossible because these `T, R, C` 192 // allows us to desambiguate a lot of configurations. 193 _phantoms: PhantomData<(T, R, C)>, 194 } 195 196 impl<T, R: Dim, C: Dim, S: fmt::Debug> fmt::Debug for Matrix<T, R, C, S> { fmt(&self, formatter: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error>197 fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> { 198 formatter 199 .debug_struct("Matrix") 200 .field("data", &self.data) 201 .finish() 202 } 203 } 204 205 impl<T, R, C, S> Default for Matrix<T, R, C, S> 206 where 207 T: Scalar, 208 R: Dim, 209 C: Dim, 210 S: Default, 211 { default() -> Self212 fn default() -> Self { 213 Matrix { 214 data: Default::default(), 215 _phantoms: PhantomData, 216 } 217 } 218 } 219 220 #[cfg(feature = "serde-serialize-no-std")] 221 impl<T, R, C, S> Serialize for Matrix<T, R, C, S> 222 where 223 T: Scalar, 224 R: Dim, 225 C: Dim, 226 S: Serialize, 227 { serialize<Ser>(&self, serializer: Ser) -> Result<Ser::Ok, Ser::Error> where Ser: Serializer,228 fn serialize<Ser>(&self, serializer: Ser) -> Result<Ser::Ok, Ser::Error> 229 where 230 Ser: Serializer, 231 { 232 self.data.serialize(serializer) 233 } 234 } 235 236 #[cfg(feature = "serde-serialize-no-std")] 237 impl<'de, T, R, C, S> Deserialize<'de> for Matrix<T, R, C, S> 238 where 239 T: Scalar, 240 R: Dim, 241 C: Dim, 242 S: Deserialize<'de>, 243 { deserialize<D>(deserializer: D) -> Result<Self, D::Error> where D: Deserializer<'de>,244 fn deserialize<D>(deserializer: D) -> Result<Self, D::Error> 245 where 246 D: Deserializer<'de>, 247 { 248 S::deserialize(deserializer).map(|x| Matrix { 249 data: x, 250 _phantoms: PhantomData, 251 }) 252 } 253 } 254 255 #[cfg(feature = "abomonation-serialize")] 256 impl<T: Scalar, R: Dim, C: Dim, S: Abomonation> Abomonation for Matrix<T, R, C, S> { entomb<W: Write>(&self, writer: &mut W) -> IOResult<()>257 unsafe fn entomb<W: Write>(&self, writer: &mut W) -> IOResult<()> { 258 self.data.entomb(writer) 259 } 260 exhume<'a, 'b>(&'a mut self, bytes: &'b mut [u8]) -> Option<&'b mut [u8]>261 unsafe fn exhume<'a, 'b>(&'a mut self, bytes: &'b mut [u8]) -> Option<&'b mut [u8]> { 262 self.data.exhume(bytes) 263 } 264 extent(&self) -> usize265 fn extent(&self) -> usize { 266 self.data.extent() 267 } 268 } 269 270 #[cfg(feature = "compare")] 271 impl<T: Scalar, R: Dim, C: Dim, S: RawStorage<T, R, C>> matrixcompare_core::Matrix<T> 272 for Matrix<T, R, C, S> 273 { rows(&self) -> usize274 fn rows(&self) -> usize { 275 self.nrows() 276 } 277 cols(&self) -> usize278 fn cols(&self) -> usize { 279 self.ncols() 280 } 281 access(&self) -> matrixcompare_core::Access<'_, T>282 fn access(&self) -> matrixcompare_core::Access<'_, T> { 283 matrixcompare_core::Access::Dense(self) 284 } 285 } 286 287 #[cfg(feature = "compare")] 288 impl<T: Scalar, R: Dim, C: Dim, S: RawStorage<T, R, C>> matrixcompare_core::DenseAccess<T> 289 for Matrix<T, R, C, S> 290 { fetch_single(&self, row: usize, col: usize) -> T291 fn fetch_single(&self, row: usize, col: usize) -> T { 292 self.index((row, col)).clone() 293 } 294 } 295 296 #[cfg(feature = "bytemuck")] 297 unsafe impl<T: Scalar, R: Dim, C: Dim, S: RawStorage<T, R, C>> bytemuck::Zeroable 298 for Matrix<T, R, C, S> 299 where 300 S: bytemuck::Zeroable, 301 { 302 } 303 304 #[cfg(feature = "bytemuck")] 305 unsafe impl<T: Scalar, R: Dim, C: Dim, S: RawStorage<T, R, C>> bytemuck::Pod for Matrix<T, R, C, S> 306 where 307 S: bytemuck::Pod, 308 Self: Copy, 309 { 310 } 311 312 #[cfg(feature = "rkyv-serialize-no-std")] 313 mod rkyv_impl { 314 use super::Matrix; 315 use core::marker::PhantomData; 316 use rkyv::{offset_of, project_struct, Archive, Deserialize, Fallible, Serialize}; 317 318 impl<T: Archive, R: Archive, C: Archive, S: Archive> Archive for Matrix<T, R, C, S> { 319 type Archived = Matrix<T::Archived, R::Archived, C::Archived, S::Archived>; 320 type Resolver = S::Resolver; 321 resolve( &self, pos: usize, resolver: Self::Resolver, out: &mut core::mem::MaybeUninit<Self::Archived>, )322 fn resolve( 323 &self, 324 pos: usize, 325 resolver: Self::Resolver, 326 out: &mut core::mem::MaybeUninit<Self::Archived>, 327 ) { 328 self.data.resolve( 329 pos + offset_of!(Self::Archived, data), 330 resolver, 331 project_struct!(out: Self::Archived => data), 332 ); 333 } 334 } 335 336 impl<T: Archive, R: Archive, C: Archive, S: Serialize<_S>, _S: Fallible + ?Sized> Serialize<_S> 337 for Matrix<T, R, C, S> 338 { serialize(&self, serializer: &mut _S) -> Result<Self::Resolver, _S::Error>339 fn serialize(&self, serializer: &mut _S) -> Result<Self::Resolver, _S::Error> { 340 self.data.serialize(serializer) 341 } 342 } 343 344 impl<T: Archive, R: Archive, C: Archive, S: Archive, D: Fallible + ?Sized> 345 Deserialize<Matrix<T, R, C, S>, D> 346 for Matrix<T::Archived, R::Archived, C::Archived, S::Archived> 347 where 348 S::Archived: Deserialize<S, D>, 349 { deserialize(&self, deserializer: &mut D) -> Result<Matrix<T, R, C, S>, D::Error>350 fn deserialize(&self, deserializer: &mut D) -> Result<Matrix<T, R, C, S>, D::Error> { 351 Ok(Matrix { 352 data: self.data.deserialize(deserializer)?, 353 _phantoms: PhantomData, 354 }) 355 } 356 } 357 } 358 359 impl<T, R, C, S> Matrix<T, R, C, S> { 360 /// Creates a new matrix with the given data without statically checking that the matrix 361 /// dimension matches the storage dimension. 362 #[inline(always)] from_data_statically_unchecked(data: S) -> Matrix<T, R, C, S>363 pub const unsafe fn from_data_statically_unchecked(data: S) -> Matrix<T, R, C, S> { 364 Matrix { 365 data, 366 _phantoms: PhantomData, 367 } 368 } 369 } 370 371 impl<T, const R: usize, const C: usize> SMatrix<T, R, C> { 372 /// Creates a new statically-allocated matrix from the given [`ArrayStorage`]. 373 /// 374 /// This method exists primarily as a workaround for the fact that `from_data` can not 375 /// work in `const fn` contexts. 376 #[inline(always)] from_array_storage(storage: ArrayStorage<T, R, C>) -> Self377 pub const fn from_array_storage(storage: ArrayStorage<T, R, C>) -> Self { 378 // This is sound because the row and column types are exactly the same as that of the 379 // storage, so there can be no mismatch 380 unsafe { Self::from_data_statically_unchecked(storage) } 381 } 382 } 383 384 // TODO: Consider removing/deprecating `from_vec_storage` once we are able to make 385 // `from_data` const fn compatible 386 #[cfg(any(feature = "std", feature = "alloc"))] 387 impl<T> DMatrix<T> { 388 /// Creates a new heap-allocated matrix from the given [`VecStorage`]. 389 /// 390 /// This method exists primarily as a workaround for the fact that `from_data` can not 391 /// work in `const fn` contexts. from_vec_storage(storage: VecStorage<T, Dynamic, Dynamic>) -> Self392 pub const fn from_vec_storage(storage: VecStorage<T, Dynamic, Dynamic>) -> Self { 393 // This is sound because the dimensions of the matrix and the storage are guaranteed 394 // to be the same 395 unsafe { Self::from_data_statically_unchecked(storage) } 396 } 397 } 398 399 // TODO: Consider removing/deprecating `from_vec_storage` once we are able to make 400 // `from_data` const fn compatible 401 #[cfg(any(feature = "std", feature = "alloc"))] 402 impl<T> DVector<T> { 403 /// Creates a new heap-allocated matrix from the given [`VecStorage`]. 404 /// 405 /// This method exists primarily as a workaround for the fact that `from_data` can not 406 /// work in `const fn` contexts. from_vec_storage(storage: VecStorage<T, Dynamic, U1>) -> Self407 pub const fn from_vec_storage(storage: VecStorage<T, Dynamic, U1>) -> Self { 408 // This is sound because the dimensions of the matrix and the storage are guaranteed 409 // to be the same 410 unsafe { Self::from_data_statically_unchecked(storage) } 411 } 412 } 413 414 impl<T, R: Dim, C: Dim> UninitMatrix<T, R, C> 415 where 416 DefaultAllocator: Allocator<T, R, C>, 417 { 418 /// Assumes a matrix's entries to be initialized. This operation should be near zero-cost. 419 /// 420 /// For the similar method that operates on matrix slices, see [`slice_assume_init`]. 421 /// 422 /// # Safety 423 /// The user must make sure that every single entry of the buffer has been initialized, 424 /// or Undefined Behavior will immediately occur. 425 #[inline(always)] assume_init(self) -> OMatrix<T, R, C>426 pub unsafe fn assume_init(self) -> OMatrix<T, R, C> { 427 OMatrix::from_data(<DefaultAllocator as Allocator<T, R, C>>::assume_init( 428 self.data, 429 )) 430 } 431 } 432 433 impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> { 434 /// Creates a new matrix with the given data. 435 #[inline(always)] from_data(data: S) -> Self436 pub fn from_data(data: S) -> Self { 437 unsafe { Self::from_data_statically_unchecked(data) } 438 } 439 440 /// The shape of this matrix returned as the tuple (number of rows, number of columns). 441 /// 442 /// # Examples: 443 /// 444 /// ``` 445 /// # use nalgebra::Matrix3x4; 446 /// let mat = Matrix3x4::<f32>::zeros(); 447 /// assert_eq!(mat.shape(), (3, 4)); 448 #[inline] 449 #[must_use] shape(&self) -> (usize, usize)450 pub fn shape(&self) -> (usize, usize) { 451 let (nrows, ncols) = self.shape_generic(); 452 (nrows.value(), ncols.value()) 453 } 454 455 /// The shape of this matrix wrapped into their representative types (`Const` or `Dynamic`). 456 #[inline] 457 #[must_use] shape_generic(&self) -> (R, C)458 pub fn shape_generic(&self) -> (R, C) { 459 self.data.shape() 460 } 461 462 /// The number of rows of this matrix. 463 /// 464 /// # Examples: 465 /// 466 /// ``` 467 /// # use nalgebra::Matrix3x4; 468 /// let mat = Matrix3x4::<f32>::zeros(); 469 /// assert_eq!(mat.nrows(), 3); 470 #[inline] 471 #[must_use] nrows(&self) -> usize472 pub fn nrows(&self) -> usize { 473 self.shape().0 474 } 475 476 /// The number of columns of this matrix. 477 /// 478 /// # Examples: 479 /// 480 /// ``` 481 /// # use nalgebra::Matrix3x4; 482 /// let mat = Matrix3x4::<f32>::zeros(); 483 /// assert_eq!(mat.ncols(), 4); 484 #[inline] 485 #[must_use] ncols(&self) -> usize486 pub fn ncols(&self) -> usize { 487 self.shape().1 488 } 489 490 /// The strides (row stride, column stride) of this matrix. 491 /// 492 /// # Examples: 493 /// 494 /// ``` 495 /// # use nalgebra::DMatrix; 496 /// let mat = DMatrix::<f32>::zeros(10, 10); 497 /// let slice = mat.slice_with_steps((0, 0), (5, 3), (1, 2)); 498 /// // The column strides is the number of steps (here 2) multiplied by the corresponding dimension. 499 /// assert_eq!(mat.strides(), (1, 10)); 500 #[inline] 501 #[must_use] strides(&self) -> (usize, usize)502 pub fn strides(&self) -> (usize, usize) { 503 let (srows, scols) = self.data.strides(); 504 (srows.value(), scols.value()) 505 } 506 507 /// Computes the row and column coordinates of the i-th element of this matrix seen as a 508 /// vector. 509 /// 510 /// # Example 511 /// ``` 512 /// # use nalgebra::Matrix2; 513 /// let m = Matrix2::new(1, 2, 514 /// 3, 4); 515 /// let i = m.vector_to_matrix_index(3); 516 /// assert_eq!(i, (1, 1)); 517 /// assert_eq!(m[i], m[3]); 518 /// ``` 519 #[inline] 520 #[must_use] vector_to_matrix_index(&self, i: usize) -> (usize, usize)521 pub fn vector_to_matrix_index(&self, i: usize) -> (usize, usize) { 522 let (nrows, ncols) = self.shape(); 523 524 // Two most common uses that should be optimized by the compiler for statically-sized 525 // matrices. 526 if nrows == 1 { 527 (0, i) 528 } else if ncols == 1 { 529 (i, 0) 530 } else { 531 (i % nrows, i / nrows) 532 } 533 } 534 535 /// Returns a pointer to the start of the matrix. 536 /// 537 /// If the matrix is not empty, this pointer is guaranteed to be aligned 538 /// and non-null. 539 /// 540 /// # Example 541 /// ``` 542 /// # use nalgebra::Matrix2; 543 /// let m = Matrix2::new(1, 2, 544 /// 3, 4); 545 /// let ptr = m.as_ptr(); 546 /// assert_eq!(unsafe { *ptr }, m[0]); 547 /// ``` 548 #[inline] 549 #[must_use] as_ptr(&self) -> *const T550 pub fn as_ptr(&self) -> *const T { 551 self.data.ptr() 552 } 553 554 /// Tests whether `self` and `rhs` are equal up to a given epsilon. 555 /// 556 /// See `relative_eq` from the `RelativeEq` trait for more details. 557 #[inline] 558 #[must_use] relative_eq<R2, C2, SB>( &self, other: &Matrix<T, R2, C2, SB>, eps: T::Epsilon, max_relative: T::Epsilon, ) -> bool where T: RelativeEq, R2: Dim, C2: Dim, SB: Storage<T, R2, C2>, T::Epsilon: Clone, ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,559 pub fn relative_eq<R2, C2, SB>( 560 &self, 561 other: &Matrix<T, R2, C2, SB>, 562 eps: T::Epsilon, 563 max_relative: T::Epsilon, 564 ) -> bool 565 where 566 T: RelativeEq, 567 R2: Dim, 568 C2: Dim, 569 SB: Storage<T, R2, C2>, 570 T::Epsilon: Clone, 571 ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>, 572 { 573 assert!(self.shape() == other.shape()); 574 self.iter() 575 .zip(other.iter()) 576 .all(|(a, b)| a.relative_eq(b, eps.clone(), max_relative.clone())) 577 } 578 579 /// Tests whether `self` and `rhs` are exactly equal. 580 #[inline] 581 #[must_use] 582 #[allow(clippy::should_implement_trait)] eq<R2, C2, SB>(&self, other: &Matrix<T, R2, C2, SB>) -> bool where T: PartialEq, R2: Dim, C2: Dim, SB: RawStorage<T, R2, C2>, ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,583 pub fn eq<R2, C2, SB>(&self, other: &Matrix<T, R2, C2, SB>) -> bool 584 where 585 T: PartialEq, 586 R2: Dim, 587 C2: Dim, 588 SB: RawStorage<T, R2, C2>, 589 ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>, 590 { 591 assert!(self.shape() == other.shape()); 592 self.iter().zip(other.iter()).all(|(a, b)| *a == *b) 593 } 594 595 /// Moves this matrix into one that owns its data. 596 #[inline] into_owned(self) -> OMatrix<T, R, C> where T: Scalar, S: Storage<T, R, C>, DefaultAllocator: Allocator<T, R, C>,597 pub fn into_owned(self) -> OMatrix<T, R, C> 598 where 599 T: Scalar, 600 S: Storage<T, R, C>, 601 DefaultAllocator: Allocator<T, R, C>, 602 { 603 Matrix::from_data(self.data.into_owned()) 604 } 605 606 // TODO: this could probably benefit from specialization. 607 // XXX: bad name. 608 /// Moves this matrix into one that owns its data. The actual type of the result depends on 609 /// matrix storage combination rules for addition. 610 #[inline] into_owned_sum<R2, C2>(self) -> MatrixSum<T, R, C, R2, C2> where T: Scalar, S: Storage<T, R, C>, R2: Dim, C2: Dim, DefaultAllocator: SameShapeAllocator<T, R, C, R2, C2>, ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,611 pub fn into_owned_sum<R2, C2>(self) -> MatrixSum<T, R, C, R2, C2> 612 where 613 T: Scalar, 614 S: Storage<T, R, C>, 615 R2: Dim, 616 C2: Dim, 617 DefaultAllocator: SameShapeAllocator<T, R, C, R2, C2>, 618 ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>, 619 { 620 if TypeId::of::<SameShapeStorage<T, R, C, R2, C2>>() == TypeId::of::<Owned<T, R, C>>() { 621 // We can just return `self.into_owned()`. 622 623 unsafe { 624 // TODO: check that those copies are optimized away by the compiler. 625 let owned = self.into_owned(); 626 let res = mem::transmute_copy(&owned); 627 mem::forget(owned); 628 res 629 } 630 } else { 631 self.clone_owned_sum() 632 } 633 } 634 635 /// Clones this matrix to one that owns its data. 636 #[inline] 637 #[must_use] clone_owned(&self) -> OMatrix<T, R, C> where T: Scalar, S: Storage<T, R, C>, DefaultAllocator: Allocator<T, R, C>,638 pub fn clone_owned(&self) -> OMatrix<T, R, C> 639 where 640 T: Scalar, 641 S: Storage<T, R, C>, 642 DefaultAllocator: Allocator<T, R, C>, 643 { 644 Matrix::from_data(self.data.clone_owned()) 645 } 646 647 /// Clones this matrix into one that owns its data. The actual type of the result depends on 648 /// matrix storage combination rules for addition. 649 #[inline] 650 #[must_use] clone_owned_sum<R2, C2>(&self) -> MatrixSum<T, R, C, R2, C2> where T: Scalar, S: Storage<T, R, C>, R2: Dim, C2: Dim, DefaultAllocator: SameShapeAllocator<T, R, C, R2, C2>, ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,651 pub fn clone_owned_sum<R2, C2>(&self) -> MatrixSum<T, R, C, R2, C2> 652 where 653 T: Scalar, 654 S: Storage<T, R, C>, 655 R2: Dim, 656 C2: Dim, 657 DefaultAllocator: SameShapeAllocator<T, R, C, R2, C2>, 658 ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>, 659 { 660 let (nrows, ncols) = self.shape(); 661 let nrows: SameShapeR<R, R2> = Dim::from_usize(nrows); 662 let ncols: SameShapeC<C, C2> = Dim::from_usize(ncols); 663 664 let mut res = Matrix::uninit(nrows, ncols); 665 666 unsafe { 667 // TODO: use copy_from? 668 for j in 0..res.ncols() { 669 for i in 0..res.nrows() { 670 *res.get_unchecked_mut((i, j)) = 671 MaybeUninit::new(self.get_unchecked((i, j)).clone()); 672 } 673 } 674 675 // SAFETY: the output has been initialized above. 676 res.assume_init() 677 } 678 } 679 680 /// Transposes `self` and store the result into `out`. 681 #[inline] transpose_to_uninit<Status, R2, C2, SB>( &self, status: Status, out: &mut Matrix<Status::Value, R2, C2, SB>, ) where Status: InitStatus<T>, T: Scalar, R2: Dim, C2: Dim, SB: RawStorageMut<Status::Value, R2, C2>, ShapeConstraint: SameNumberOfRows<R, C2> + SameNumberOfColumns<C, R2>,682 fn transpose_to_uninit<Status, R2, C2, SB>( 683 &self, 684 status: Status, 685 out: &mut Matrix<Status::Value, R2, C2, SB>, 686 ) where 687 Status: InitStatus<T>, 688 T: Scalar, 689 R2: Dim, 690 C2: Dim, 691 SB: RawStorageMut<Status::Value, R2, C2>, 692 ShapeConstraint: SameNumberOfRows<R, C2> + SameNumberOfColumns<C, R2>, 693 { 694 let (nrows, ncols) = self.shape(); 695 assert!( 696 (ncols, nrows) == out.shape(), 697 "Incompatible shape for transposition." 698 ); 699 700 // TODO: optimize that. 701 for i in 0..nrows { 702 for j in 0..ncols { 703 // Safety: the indices are in range. 704 unsafe { 705 Status::init( 706 out.get_unchecked_mut((j, i)), 707 self.get_unchecked((i, j)).clone(), 708 ); 709 } 710 } 711 } 712 } 713 714 /// Transposes `self` and store the result into `out`. 715 #[inline] transpose_to<R2, C2, SB>(&self, out: &mut Matrix<T, R2, C2, SB>) where T: Scalar, R2: Dim, C2: Dim, SB: RawStorageMut<T, R2, C2>, ShapeConstraint: SameNumberOfRows<R, C2> + SameNumberOfColumns<C, R2>,716 pub fn transpose_to<R2, C2, SB>(&self, out: &mut Matrix<T, R2, C2, SB>) 717 where 718 T: Scalar, 719 R2: Dim, 720 C2: Dim, 721 SB: RawStorageMut<T, R2, C2>, 722 ShapeConstraint: SameNumberOfRows<R, C2> + SameNumberOfColumns<C, R2>, 723 { 724 self.transpose_to_uninit(Init, out) 725 } 726 727 /// Transposes `self`. 728 #[inline] 729 #[must_use = "Did you mean to use transpose_mut()?"] transpose(&self) -> OMatrix<T, C, R> where T: Scalar, DefaultAllocator: Allocator<T, C, R>,730 pub fn transpose(&self) -> OMatrix<T, C, R> 731 where 732 T: Scalar, 733 DefaultAllocator: Allocator<T, C, R>, 734 { 735 let (nrows, ncols) = self.shape_generic(); 736 737 let mut res = Matrix::uninit(ncols, nrows); 738 self.transpose_to_uninit(Uninit, &mut res); 739 // Safety: res is now fully initialized. 740 unsafe { res.assume_init() } 741 } 742 } 743 744 /// # Elementwise mapping and folding 745 impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> { 746 /// Returns a matrix containing the result of `f` applied to each of its entries. 747 #[inline] 748 #[must_use] map<T2: Scalar, F: FnMut(T) -> T2>(&self, mut f: F) -> OMatrix<T2, R, C> where T: Scalar, DefaultAllocator: Allocator<T2, R, C>,749 pub fn map<T2: Scalar, F: FnMut(T) -> T2>(&self, mut f: F) -> OMatrix<T2, R, C> 750 where 751 T: Scalar, 752 DefaultAllocator: Allocator<T2, R, C>, 753 { 754 let (nrows, ncols) = self.shape_generic(); 755 let mut res = Matrix::uninit(nrows, ncols); 756 757 for j in 0..ncols.value() { 758 for i in 0..nrows.value() { 759 // Safety: all indices are in range. 760 unsafe { 761 let a = self.data.get_unchecked(i, j).clone(); 762 *res.data.get_unchecked_mut(i, j) = MaybeUninit::new(f(a)); 763 } 764 } 765 } 766 767 // Safety: res is now fully initialized. 768 unsafe { res.assume_init() } 769 } 770 771 /// Cast the components of `self` to another type. 772 /// 773 /// # Example 774 /// ``` 775 /// # use nalgebra::Vector3; 776 /// let q = Vector3::new(1.0f64, 2.0, 3.0); 777 /// let q2 = q.cast::<f32>(); 778 /// assert_eq!(q2, Vector3::new(1.0f32, 2.0, 3.0)); 779 /// ``` cast<T2: Scalar>(self) -> OMatrix<T2, R, C> where T: Scalar, OMatrix<T2, R, C>: SupersetOf<Self>, DefaultAllocator: Allocator<T2, R, C>,780 pub fn cast<T2: Scalar>(self) -> OMatrix<T2, R, C> 781 where 782 T: Scalar, 783 OMatrix<T2, R, C>: SupersetOf<Self>, 784 DefaultAllocator: Allocator<T2, R, C>, 785 { 786 crate::convert(self) 787 } 788 789 /// Similar to `self.iter().fold(init, f)` except that `init` is replaced by a closure. 790 /// 791 /// The initialization closure is given the first component of this matrix: 792 /// - If the matrix has no component (0 rows or 0 columns) then `init_f` is called with `None` 793 /// and its return value is the value returned by this method. 794 /// - If the matrix has has least one component, then `init_f` is called with the first component 795 /// to compute the initial value. Folding then continues on all the remaining components of the matrix. 796 #[inline] 797 #[must_use] fold_with<T2>( &self, init_f: impl FnOnce(Option<&T>) -> T2, f: impl FnMut(T2, &T) -> T2, ) -> T2 where T: Scalar,798 pub fn fold_with<T2>( 799 &self, 800 init_f: impl FnOnce(Option<&T>) -> T2, 801 f: impl FnMut(T2, &T) -> T2, 802 ) -> T2 803 where 804 T: Scalar, 805 { 806 let mut it = self.iter(); 807 let init = init_f(it.next()); 808 it.fold(init, f) 809 } 810 811 /// Returns a matrix containing the result of `f` applied to each of its entries. Unlike `map`, 812 /// `f` also gets passed the row and column index, i.e. `f(row, col, value)`. 813 #[inline] 814 #[must_use] map_with_location<T2: Scalar, F: FnMut(usize, usize, T) -> T2>( &self, mut f: F, ) -> OMatrix<T2, R, C> where T: Scalar, DefaultAllocator: Allocator<T2, R, C>,815 pub fn map_with_location<T2: Scalar, F: FnMut(usize, usize, T) -> T2>( 816 &self, 817 mut f: F, 818 ) -> OMatrix<T2, R, C> 819 where 820 T: Scalar, 821 DefaultAllocator: Allocator<T2, R, C>, 822 { 823 let (nrows, ncols) = self.shape_generic(); 824 let mut res = Matrix::uninit(nrows, ncols); 825 826 for j in 0..ncols.value() { 827 for i in 0..nrows.value() { 828 // Safety: all indices are in range. 829 unsafe { 830 let a = self.data.get_unchecked(i, j).clone(); 831 *res.data.get_unchecked_mut(i, j) = MaybeUninit::new(f(i, j, a)); 832 } 833 } 834 } 835 836 // Safety: res is now fully initialized. 837 unsafe { res.assume_init() } 838 } 839 840 /// Returns a matrix containing the result of `f` applied to each entries of `self` and 841 /// `rhs`. 842 #[inline] 843 #[must_use] zip_map<T2, N3, S2, F>(&self, rhs: &Matrix<T2, R, C, S2>, mut f: F) -> OMatrix<N3, R, C> where T: Scalar, T2: Scalar, N3: Scalar, S2: RawStorage<T2, R, C>, F: FnMut(T, T2) -> N3, DefaultAllocator: Allocator<N3, R, C>,844 pub fn zip_map<T2, N3, S2, F>(&self, rhs: &Matrix<T2, R, C, S2>, mut f: F) -> OMatrix<N3, R, C> 845 where 846 T: Scalar, 847 T2: Scalar, 848 N3: Scalar, 849 S2: RawStorage<T2, R, C>, 850 F: FnMut(T, T2) -> N3, 851 DefaultAllocator: Allocator<N3, R, C>, 852 { 853 let (nrows, ncols) = self.shape_generic(); 854 let mut res = Matrix::uninit(nrows, ncols); 855 856 assert_eq!( 857 (nrows.value(), ncols.value()), 858 rhs.shape(), 859 "Matrix simultaneous traversal error: dimension mismatch." 860 ); 861 862 for j in 0..ncols.value() { 863 for i in 0..nrows.value() { 864 // Safety: all indices are in range. 865 unsafe { 866 let a = self.data.get_unchecked(i, j).clone(); 867 let b = rhs.data.get_unchecked(i, j).clone(); 868 *res.data.get_unchecked_mut(i, j) = MaybeUninit::new(f(a, b)) 869 } 870 } 871 } 872 873 // Safety: res is now fully initialized. 874 unsafe { res.assume_init() } 875 } 876 877 /// Returns a matrix containing the result of `f` applied to each entries of `self` and 878 /// `b`, and `c`. 879 #[inline] 880 #[must_use] zip_zip_map<T2, N3, N4, S2, S3, F>( &self, b: &Matrix<T2, R, C, S2>, c: &Matrix<N3, R, C, S3>, mut f: F, ) -> OMatrix<N4, R, C> where T: Scalar, T2: Scalar, N3: Scalar, N4: Scalar, S2: RawStorage<T2, R, C>, S3: RawStorage<N3, R, C>, F: FnMut(T, T2, N3) -> N4, DefaultAllocator: Allocator<N4, R, C>,881 pub fn zip_zip_map<T2, N3, N4, S2, S3, F>( 882 &self, 883 b: &Matrix<T2, R, C, S2>, 884 c: &Matrix<N3, R, C, S3>, 885 mut f: F, 886 ) -> OMatrix<N4, R, C> 887 where 888 T: Scalar, 889 T2: Scalar, 890 N3: Scalar, 891 N4: Scalar, 892 S2: RawStorage<T2, R, C>, 893 S3: RawStorage<N3, R, C>, 894 F: FnMut(T, T2, N3) -> N4, 895 DefaultAllocator: Allocator<N4, R, C>, 896 { 897 let (nrows, ncols) = self.shape_generic(); 898 let mut res = Matrix::uninit(nrows, ncols); 899 900 assert_eq!( 901 (nrows.value(), ncols.value()), 902 b.shape(), 903 "Matrix simultaneous traversal error: dimension mismatch." 904 ); 905 assert_eq!( 906 (nrows.value(), ncols.value()), 907 c.shape(), 908 "Matrix simultaneous traversal error: dimension mismatch." 909 ); 910 911 for j in 0..ncols.value() { 912 for i in 0..nrows.value() { 913 // Safety: all indices are in range. 914 unsafe { 915 let a = self.data.get_unchecked(i, j).clone(); 916 let b = b.data.get_unchecked(i, j).clone(); 917 let c = c.data.get_unchecked(i, j).clone(); 918 *res.data.get_unchecked_mut(i, j) = MaybeUninit::new(f(a, b, c)) 919 } 920 } 921 } 922 923 // Safety: res is now fully initialized. 924 unsafe { res.assume_init() } 925 } 926 927 /// Folds a function `f` on each entry of `self`. 928 #[inline] 929 #[must_use] fold<Acc>(&self, init: Acc, mut f: impl FnMut(Acc, T) -> Acc) -> Acc where T: Scalar,930 pub fn fold<Acc>(&self, init: Acc, mut f: impl FnMut(Acc, T) -> Acc) -> Acc 931 where 932 T: Scalar, 933 { 934 let (nrows, ncols) = self.shape_generic(); 935 936 let mut res = init; 937 938 for j in 0..ncols.value() { 939 for i in 0..nrows.value() { 940 // Safety: all indices are in range. 941 unsafe { 942 let a = self.data.get_unchecked(i, j).clone(); 943 res = f(res, a) 944 } 945 } 946 } 947 948 res 949 } 950 951 /// Folds a function `f` on each pairs of entries from `self` and `rhs`. 952 #[inline] 953 #[must_use] zip_fold<T2, R2, C2, S2, Acc>( &self, rhs: &Matrix<T2, R2, C2, S2>, init: Acc, mut f: impl FnMut(Acc, T, T2) -> Acc, ) -> Acc where T: Scalar, T2: Scalar, R2: Dim, C2: Dim, S2: RawStorage<T2, R2, C2>, ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,954 pub fn zip_fold<T2, R2, C2, S2, Acc>( 955 &self, 956 rhs: &Matrix<T2, R2, C2, S2>, 957 init: Acc, 958 mut f: impl FnMut(Acc, T, T2) -> Acc, 959 ) -> Acc 960 where 961 T: Scalar, 962 T2: Scalar, 963 R2: Dim, 964 C2: Dim, 965 S2: RawStorage<T2, R2, C2>, 966 ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>, 967 { 968 let (nrows, ncols) = self.shape_generic(); 969 970 let mut res = init; 971 972 assert_eq!( 973 (nrows.value(), ncols.value()), 974 rhs.shape(), 975 "Matrix simultaneous traversal error: dimension mismatch." 976 ); 977 978 for j in 0..ncols.value() { 979 for i in 0..nrows.value() { 980 unsafe { 981 let a = self.data.get_unchecked(i, j).clone(); 982 let b = rhs.data.get_unchecked(i, j).clone(); 983 res = f(res, a, b) 984 } 985 } 986 } 987 988 res 989 } 990 991 /// Applies a closure `f` to modify each component of `self`. 992 #[inline] apply<F: FnMut(&mut T)>(&mut self, mut f: F) where S: RawStorageMut<T, R, C>,993 pub fn apply<F: FnMut(&mut T)>(&mut self, mut f: F) 994 where 995 S: RawStorageMut<T, R, C>, 996 { 997 let (nrows, ncols) = self.shape(); 998 999 for j in 0..ncols { 1000 for i in 0..nrows { 1001 unsafe { 1002 let e = self.data.get_unchecked_mut(i, j); 1003 f(e) 1004 } 1005 } 1006 } 1007 } 1008 1009 /// Replaces each component of `self` by the result of a closure `f` applied on its components 1010 /// joined with the components from `rhs`. 1011 #[inline] zip_apply<T2, R2, C2, S2>( &mut self, rhs: &Matrix<T2, R2, C2, S2>, mut f: impl FnMut(&mut T, T2), ) where S: RawStorageMut<T, R, C>, T2: Scalar, R2: Dim, C2: Dim, S2: RawStorage<T2, R2, C2>, ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,1012 pub fn zip_apply<T2, R2, C2, S2>( 1013 &mut self, 1014 rhs: &Matrix<T2, R2, C2, S2>, 1015 mut f: impl FnMut(&mut T, T2), 1016 ) where 1017 S: RawStorageMut<T, R, C>, 1018 T2: Scalar, 1019 R2: Dim, 1020 C2: Dim, 1021 S2: RawStorage<T2, R2, C2>, 1022 ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>, 1023 { 1024 let (nrows, ncols) = self.shape(); 1025 1026 assert_eq!( 1027 (nrows, ncols), 1028 rhs.shape(), 1029 "Matrix simultaneous traversal error: dimension mismatch." 1030 ); 1031 1032 for j in 0..ncols { 1033 for i in 0..nrows { 1034 unsafe { 1035 let e = self.data.get_unchecked_mut(i, j); 1036 let rhs = rhs.get_unchecked((i, j)).clone(); 1037 f(e, rhs) 1038 } 1039 } 1040 } 1041 } 1042 1043 /// Replaces each component of `self` by the result of a closure `f` applied on its components 1044 /// joined with the components from `b` and `c`. 1045 #[inline] zip_zip_apply<T2, R2, C2, S2, N3, R3, C3, S3>( &mut self, b: &Matrix<T2, R2, C2, S2>, c: &Matrix<N3, R3, C3, S3>, mut f: impl FnMut(&mut T, T2, N3), ) where S: RawStorageMut<T, R, C>, T2: Scalar, R2: Dim, C2: Dim, S2: RawStorage<T2, R2, C2>, N3: Scalar, R3: Dim, C3: Dim, S3: RawStorage<N3, R3, C3>, ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>, ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,1046 pub fn zip_zip_apply<T2, R2, C2, S2, N3, R3, C3, S3>( 1047 &mut self, 1048 b: &Matrix<T2, R2, C2, S2>, 1049 c: &Matrix<N3, R3, C3, S3>, 1050 mut f: impl FnMut(&mut T, T2, N3), 1051 ) where 1052 S: RawStorageMut<T, R, C>, 1053 T2: Scalar, 1054 R2: Dim, 1055 C2: Dim, 1056 S2: RawStorage<T2, R2, C2>, 1057 N3: Scalar, 1058 R3: Dim, 1059 C3: Dim, 1060 S3: RawStorage<N3, R3, C3>, 1061 ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>, 1062 ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>, 1063 { 1064 let (nrows, ncols) = self.shape(); 1065 1066 assert_eq!( 1067 (nrows, ncols), 1068 b.shape(), 1069 "Matrix simultaneous traversal error: dimension mismatch." 1070 ); 1071 assert_eq!( 1072 (nrows, ncols), 1073 c.shape(), 1074 "Matrix simultaneous traversal error: dimension mismatch." 1075 ); 1076 1077 for j in 0..ncols { 1078 for i in 0..nrows { 1079 unsafe { 1080 let e = self.data.get_unchecked_mut(i, j); 1081 let b = b.get_unchecked((i, j)).clone(); 1082 let c = c.get_unchecked((i, j)).clone(); 1083 f(e, b, c) 1084 } 1085 } 1086 } 1087 } 1088 } 1089 1090 /// # Iteration on components, rows, and columns 1091 impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> { 1092 /// Iterates through this matrix coordinates in column-major order. 1093 /// 1094 /// # Examples: 1095 /// 1096 /// ``` 1097 /// # use nalgebra::Matrix2x3; 1098 /// let mat = Matrix2x3::new(11, 12, 13, 1099 /// 21, 22, 23); 1100 /// let mut it = mat.iter(); 1101 /// assert_eq!(*it.next().unwrap(), 11); 1102 /// assert_eq!(*it.next().unwrap(), 21); 1103 /// assert_eq!(*it.next().unwrap(), 12); 1104 /// assert_eq!(*it.next().unwrap(), 22); 1105 /// assert_eq!(*it.next().unwrap(), 13); 1106 /// assert_eq!(*it.next().unwrap(), 23); 1107 /// assert!(it.next().is_none()); 1108 #[inline] iter(&self) -> MatrixIter<'_, T, R, C, S>1109 pub fn iter(&self) -> MatrixIter<'_, T, R, C, S> { 1110 MatrixIter::new(&self.data) 1111 } 1112 1113 /// Iterate through the rows of this matrix. 1114 /// 1115 /// # Example 1116 /// ``` 1117 /// # use nalgebra::Matrix2x3; 1118 /// let mut a = Matrix2x3::new(1, 2, 3, 1119 /// 4, 5, 6); 1120 /// for (i, row) in a.row_iter().enumerate() { 1121 /// assert_eq!(row, a.row(i)) 1122 /// } 1123 /// ``` 1124 #[inline] row_iter(&self) -> RowIter<'_, T, R, C, S>1125 pub fn row_iter(&self) -> RowIter<'_, T, R, C, S> { 1126 RowIter::new(self) 1127 } 1128 1129 /// Iterate through the columns of this matrix. 1130 /// # Example 1131 /// ``` 1132 /// # use nalgebra::Matrix2x3; 1133 /// let mut a = Matrix2x3::new(1, 2, 3, 1134 /// 4, 5, 6); 1135 /// for (i, column) in a.column_iter().enumerate() { 1136 /// assert_eq!(column, a.column(i)) 1137 /// } 1138 /// ``` 1139 #[inline] column_iter(&self) -> ColumnIter<'_, T, R, C, S>1140 pub fn column_iter(&self) -> ColumnIter<'_, T, R, C, S> { 1141 ColumnIter::new(self) 1142 } 1143 1144 /// Mutably iterates through this matrix coordinates. 1145 #[inline] iter_mut(&mut self) -> MatrixIterMut<'_, T, R, C, S> where S: RawStorageMut<T, R, C>,1146 pub fn iter_mut(&mut self) -> MatrixIterMut<'_, T, R, C, S> 1147 where 1148 S: RawStorageMut<T, R, C>, 1149 { 1150 MatrixIterMut::new(&mut self.data) 1151 } 1152 1153 /// Mutably iterates through this matrix rows. 1154 /// 1155 /// # Example 1156 /// ``` 1157 /// # use nalgebra::Matrix2x3; 1158 /// let mut a = Matrix2x3::new(1, 2, 3, 1159 /// 4, 5, 6); 1160 /// for (i, mut row) in a.row_iter_mut().enumerate() { 1161 /// row *= (i + 1) * 10; 1162 /// } 1163 /// 1164 /// let expected = Matrix2x3::new(10, 20, 30, 1165 /// 80, 100, 120); 1166 /// assert_eq!(a, expected); 1167 /// ``` 1168 #[inline] row_iter_mut(&mut self) -> RowIterMut<'_, T, R, C, S> where S: RawStorageMut<T, R, C>,1169 pub fn row_iter_mut(&mut self) -> RowIterMut<'_, T, R, C, S> 1170 where 1171 S: RawStorageMut<T, R, C>, 1172 { 1173 RowIterMut::new(self) 1174 } 1175 1176 /// Mutably iterates through this matrix columns. 1177 /// 1178 /// # Example 1179 /// ``` 1180 /// # use nalgebra::Matrix2x3; 1181 /// let mut a = Matrix2x3::new(1, 2, 3, 1182 /// 4, 5, 6); 1183 /// for (i, mut col) in a.column_iter_mut().enumerate() { 1184 /// col *= (i + 1) * 10; 1185 /// } 1186 /// 1187 /// let expected = Matrix2x3::new(10, 40, 90, 1188 /// 40, 100, 180); 1189 /// assert_eq!(a, expected); 1190 /// ``` 1191 #[inline] column_iter_mut(&mut self) -> ColumnIterMut<'_, T, R, C, S> where S: RawStorageMut<T, R, C>,1192 pub fn column_iter_mut(&mut self) -> ColumnIterMut<'_, T, R, C, S> 1193 where 1194 S: RawStorageMut<T, R, C>, 1195 { 1196 ColumnIterMut::new(self) 1197 } 1198 } 1199 1200 impl<T, R: Dim, C: Dim, S: RawStorageMut<T, R, C>> Matrix<T, R, C, S> { 1201 /// Returns a mutable pointer to the start of the matrix. 1202 /// 1203 /// If the matrix is not empty, this pointer is guaranteed to be aligned 1204 /// and non-null. 1205 #[inline] as_mut_ptr(&mut self) -> *mut T1206 pub fn as_mut_ptr(&mut self) -> *mut T { 1207 self.data.ptr_mut() 1208 } 1209 1210 /// Swaps two entries without bound-checking. 1211 #[inline] swap_unchecked(&mut self, row_cols1: (usize, usize), row_cols2: (usize, usize))1212 pub unsafe fn swap_unchecked(&mut self, row_cols1: (usize, usize), row_cols2: (usize, usize)) { 1213 debug_assert!(row_cols1.0 < self.nrows() && row_cols1.1 < self.ncols()); 1214 debug_assert!(row_cols2.0 < self.nrows() && row_cols2.1 < self.ncols()); 1215 self.data.swap_unchecked(row_cols1, row_cols2) 1216 } 1217 1218 /// Swaps two entries. 1219 #[inline] swap(&mut self, row_cols1: (usize, usize), row_cols2: (usize, usize))1220 pub fn swap(&mut self, row_cols1: (usize, usize), row_cols2: (usize, usize)) { 1221 let (nrows, ncols) = self.shape(); 1222 assert!( 1223 row_cols1.0 < nrows && row_cols1.1 < ncols, 1224 "Matrix elements swap index out of bounds." 1225 ); 1226 assert!( 1227 row_cols2.0 < nrows && row_cols2.1 < ncols, 1228 "Matrix elements swap index out of bounds." 1229 ); 1230 unsafe { self.swap_unchecked(row_cols1, row_cols2) } 1231 } 1232 1233 /// Fills this matrix with the content of a slice. Both must hold the same number of elements. 1234 /// 1235 /// The components of the slice are assumed to be ordered in column-major order. 1236 #[inline] copy_from_slice(&mut self, slice: &[T]) where T: Scalar,1237 pub fn copy_from_slice(&mut self, slice: &[T]) 1238 where 1239 T: Scalar, 1240 { 1241 let (nrows, ncols) = self.shape(); 1242 1243 assert!( 1244 nrows * ncols == slice.len(), 1245 "The slice must contain the same number of elements as the matrix." 1246 ); 1247 1248 for j in 0..ncols { 1249 for i in 0..nrows { 1250 unsafe { 1251 *self.get_unchecked_mut((i, j)) = slice.get_unchecked(i + j * nrows).clone(); 1252 } 1253 } 1254 } 1255 } 1256 1257 /// Fills this matrix with the content of another one. Both must have the same shape. 1258 #[inline] copy_from<R2, C2, SB>(&mut self, other: &Matrix<T, R2, C2, SB>) where T: Scalar, R2: Dim, C2: Dim, SB: RawStorage<T, R2, C2>, ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,1259 pub fn copy_from<R2, C2, SB>(&mut self, other: &Matrix<T, R2, C2, SB>) 1260 where 1261 T: Scalar, 1262 R2: Dim, 1263 C2: Dim, 1264 SB: RawStorage<T, R2, C2>, 1265 ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>, 1266 { 1267 assert!( 1268 self.shape() == other.shape(), 1269 "Unable to copy from a matrix with a different shape." 1270 ); 1271 1272 for j in 0..self.ncols() { 1273 for i in 0..self.nrows() { 1274 unsafe { 1275 *self.get_unchecked_mut((i, j)) = other.get_unchecked((i, j)).clone(); 1276 } 1277 } 1278 } 1279 } 1280 1281 /// Fills this matrix with the content of the transpose another one. 1282 #[inline] tr_copy_from<R2, C2, SB>(&mut self, other: &Matrix<T, R2, C2, SB>) where T: Scalar, R2: Dim, C2: Dim, SB: RawStorage<T, R2, C2>, ShapeConstraint: DimEq<R, C2> + SameNumberOfColumns<C, R2>,1283 pub fn tr_copy_from<R2, C2, SB>(&mut self, other: &Matrix<T, R2, C2, SB>) 1284 where 1285 T: Scalar, 1286 R2: Dim, 1287 C2: Dim, 1288 SB: RawStorage<T, R2, C2>, 1289 ShapeConstraint: DimEq<R, C2> + SameNumberOfColumns<C, R2>, 1290 { 1291 let (nrows, ncols) = self.shape(); 1292 assert!( 1293 (ncols, nrows) == other.shape(), 1294 "Unable to copy from a matrix with incompatible shape." 1295 ); 1296 1297 for j in 0..ncols { 1298 for i in 0..nrows { 1299 unsafe { 1300 *self.get_unchecked_mut((i, j)) = other.get_unchecked((j, i)).clone(); 1301 } 1302 } 1303 } 1304 } 1305 1306 // TODO: rename `apply` to `apply_mut` and `apply_into` to `apply`? 1307 /// Returns `self` with each of its components replaced by the result of a closure `f` applied on it. 1308 #[inline] apply_into<F: FnMut(&mut T)>(mut self, f: F) -> Self1309 pub fn apply_into<F: FnMut(&mut T)>(mut self, f: F) -> Self { 1310 self.apply(f); 1311 self 1312 } 1313 } 1314 1315 impl<T, D: Dim, S: RawStorage<T, D>> Vector<T, D, S> { 1316 /// Gets a reference to the i-th element of this column vector without bound checking. 1317 #[inline] 1318 #[must_use] vget_unchecked(&self, i: usize) -> &T1319 pub unsafe fn vget_unchecked(&self, i: usize) -> &T { 1320 debug_assert!(i < self.nrows(), "Vector index out of bounds."); 1321 let i = i * self.strides().0; 1322 self.data.get_unchecked_linear(i) 1323 } 1324 } 1325 1326 impl<T, D: Dim, S: RawStorageMut<T, D>> Vector<T, D, S> { 1327 /// Gets a mutable reference to the i-th element of this column vector without bound checking. 1328 #[inline] 1329 #[must_use] vget_unchecked_mut(&mut self, i: usize) -> &mut T1330 pub unsafe fn vget_unchecked_mut(&mut self, i: usize) -> &mut T { 1331 debug_assert!(i < self.nrows(), "Vector index out of bounds."); 1332 let i = i * self.strides().0; 1333 self.data.get_unchecked_linear_mut(i) 1334 } 1335 } 1336 1337 impl<T, R: Dim, C: Dim, S: RawStorage<T, R, C> + IsContiguous> Matrix<T, R, C, S> { 1338 /// Extracts a slice containing the entire matrix entries ordered column-by-columns. 1339 #[inline] 1340 #[must_use] as_slice(&self) -> &[T]1341 pub fn as_slice(&self) -> &[T] { 1342 // Safety: this is OK thanks to the IsContiguous trait. 1343 unsafe { self.data.as_slice_unchecked() } 1344 } 1345 } 1346 1347 impl<T, R: Dim, C: Dim, S: RawStorageMut<T, R, C> + IsContiguous> Matrix<T, R, C, S> { 1348 /// Extracts a mutable slice containing the entire matrix entries ordered column-by-columns. 1349 #[inline] 1350 #[must_use] as_mut_slice(&mut self) -> &mut [T]1351 pub fn as_mut_slice(&mut self) -> &mut [T] { 1352 // Safety: this is OK thanks to the IsContiguous trait. 1353 unsafe { self.data.as_mut_slice_unchecked() } 1354 } 1355 } 1356 1357 impl<T: Scalar, D: Dim, S: RawStorageMut<T, D, D>> Matrix<T, D, D, S> { 1358 /// Transposes the square matrix `self` in-place. transpose_mut(&mut self)1359 pub fn transpose_mut(&mut self) { 1360 assert!( 1361 self.is_square(), 1362 "Unable to transpose a non-square matrix in-place." 1363 ); 1364 1365 let dim = self.shape().0; 1366 1367 for i in 1..dim { 1368 for j in 0..i { 1369 unsafe { self.swap_unchecked((i, j), (j, i)) } 1370 } 1371 } 1372 } 1373 } 1374 1375 impl<T: SimdComplexField, R: Dim, C: Dim, S: RawStorage<T, R, C>> Matrix<T, R, C, S> { 1376 /// Takes the adjoint (aka. conjugate-transpose) of `self` and store the result into `out`. 1377 #[inline] adjoint_to_uninit<Status, R2, C2, SB>( &self, status: Status, out: &mut Matrix<Status::Value, R2, C2, SB>, ) where Status: InitStatus<T>, R2: Dim, C2: Dim, SB: RawStorageMut<Status::Value, R2, C2>, ShapeConstraint: SameNumberOfRows<R, C2> + SameNumberOfColumns<C, R2>,1378 fn adjoint_to_uninit<Status, R2, C2, SB>( 1379 &self, 1380 status: Status, 1381 out: &mut Matrix<Status::Value, R2, C2, SB>, 1382 ) where 1383 Status: InitStatus<T>, 1384 R2: Dim, 1385 C2: Dim, 1386 SB: RawStorageMut<Status::Value, R2, C2>, 1387 ShapeConstraint: SameNumberOfRows<R, C2> + SameNumberOfColumns<C, R2>, 1388 { 1389 let (nrows, ncols) = self.shape(); 1390 assert!( 1391 (ncols, nrows) == out.shape(), 1392 "Incompatible shape for transpose-copy." 1393 ); 1394 1395 // TODO: optimize that. 1396 for i in 0..nrows { 1397 for j in 0..ncols { 1398 // Safety: all indices are in range. 1399 unsafe { 1400 Status::init( 1401 out.get_unchecked_mut((j, i)), 1402 self.get_unchecked((i, j)).clone().simd_conjugate(), 1403 ); 1404 } 1405 } 1406 } 1407 } 1408 1409 /// Takes the adjoint (aka. conjugate-transpose) of `self` and store the result into `out`. 1410 #[inline] adjoint_to<R2, C2, SB>(&self, out: &mut Matrix<T, R2, C2, SB>) where R2: Dim, C2: Dim, SB: RawStorageMut<T, R2, C2>, ShapeConstraint: SameNumberOfRows<R, C2> + SameNumberOfColumns<C, R2>,1411 pub fn adjoint_to<R2, C2, SB>(&self, out: &mut Matrix<T, R2, C2, SB>) 1412 where 1413 R2: Dim, 1414 C2: Dim, 1415 SB: RawStorageMut<T, R2, C2>, 1416 ShapeConstraint: SameNumberOfRows<R, C2> + SameNumberOfColumns<C, R2>, 1417 { 1418 self.adjoint_to_uninit(Init, out) 1419 } 1420 1421 /// The adjoint (aka. conjugate-transpose) of `self`. 1422 #[inline] 1423 #[must_use = "Did you mean to use adjoint_mut()?"] adjoint(&self) -> OMatrix<T, C, R> where DefaultAllocator: Allocator<T, C, R>,1424 pub fn adjoint(&self) -> OMatrix<T, C, R> 1425 where 1426 DefaultAllocator: Allocator<T, C, R>, 1427 { 1428 let (nrows, ncols) = self.shape_generic(); 1429 1430 let mut res = Matrix::uninit(ncols, nrows); 1431 self.adjoint_to_uninit(Uninit, &mut res); 1432 1433 // Safety: res is now fully initialized. 1434 unsafe { res.assume_init() } 1435 } 1436 1437 /// Takes the conjugate and transposes `self` and store the result into `out`. 1438 #[deprecated(note = "Renamed `self.adjoint_to(out)`.")] 1439 #[inline] conjugate_transpose_to<R2, C2, SB>(&self, out: &mut Matrix<T, R2, C2, SB>) where R2: Dim, C2: Dim, SB: RawStorageMut<T, R2, C2>, ShapeConstraint: SameNumberOfRows<R, C2> + SameNumberOfColumns<C, R2>,1440 pub fn conjugate_transpose_to<R2, C2, SB>(&self, out: &mut Matrix<T, R2, C2, SB>) 1441 where 1442 R2: Dim, 1443 C2: Dim, 1444 SB: RawStorageMut<T, R2, C2>, 1445 ShapeConstraint: SameNumberOfRows<R, C2> + SameNumberOfColumns<C, R2>, 1446 { 1447 self.adjoint_to(out) 1448 } 1449 1450 /// The conjugate transposition of `self`. 1451 #[deprecated(note = "Renamed `self.adjoint()`.")] 1452 #[inline] conjugate_transpose(&self) -> OMatrix<T, C, R> where DefaultAllocator: Allocator<T, C, R>,1453 pub fn conjugate_transpose(&self) -> OMatrix<T, C, R> 1454 where 1455 DefaultAllocator: Allocator<T, C, R>, 1456 { 1457 self.adjoint() 1458 } 1459 1460 /// The conjugate of `self`. 1461 #[inline] 1462 #[must_use = "Did you mean to use conjugate_mut()?"] conjugate(&self) -> OMatrix<T, R, C> where DefaultAllocator: Allocator<T, R, C>,1463 pub fn conjugate(&self) -> OMatrix<T, R, C> 1464 where 1465 DefaultAllocator: Allocator<T, R, C>, 1466 { 1467 self.map(|e| e.simd_conjugate()) 1468 } 1469 1470 /// Divides each component of the complex matrix `self` by the given real. 1471 #[inline] 1472 #[must_use = "Did you mean to use unscale_mut()?"] unscale(&self, real: T::SimdRealField) -> OMatrix<T, R, C> where DefaultAllocator: Allocator<T, R, C>,1473 pub fn unscale(&self, real: T::SimdRealField) -> OMatrix<T, R, C> 1474 where 1475 DefaultAllocator: Allocator<T, R, C>, 1476 { 1477 self.map(|e| e.simd_unscale(real.clone())) 1478 } 1479 1480 /// Multiplies each component of the complex matrix `self` by the given real. 1481 #[inline] 1482 #[must_use = "Did you mean to use scale_mut()?"] scale(&self, real: T::SimdRealField) -> OMatrix<T, R, C> where DefaultAllocator: Allocator<T, R, C>,1483 pub fn scale(&self, real: T::SimdRealField) -> OMatrix<T, R, C> 1484 where 1485 DefaultAllocator: Allocator<T, R, C>, 1486 { 1487 self.map(|e| e.simd_scale(real.clone())) 1488 } 1489 } 1490 1491 impl<T: SimdComplexField, R: Dim, C: Dim, S: RawStorageMut<T, R, C>> Matrix<T, R, C, S> { 1492 /// The conjugate of the complex matrix `self` computed in-place. 1493 #[inline] conjugate_mut(&mut self)1494 pub fn conjugate_mut(&mut self) { 1495 self.apply(|e| *e = e.clone().simd_conjugate()) 1496 } 1497 1498 /// Divides each component of the complex matrix `self` by the given real. 1499 #[inline] unscale_mut(&mut self, real: T::SimdRealField)1500 pub fn unscale_mut(&mut self, real: T::SimdRealField) { 1501 self.apply(|e| *e = e.clone().simd_unscale(real.clone())) 1502 } 1503 1504 /// Multiplies each component of the complex matrix `self` by the given real. 1505 #[inline] scale_mut(&mut self, real: T::SimdRealField)1506 pub fn scale_mut(&mut self, real: T::SimdRealField) { 1507 self.apply(|e| *e = e.clone().simd_scale(real.clone())) 1508 } 1509 } 1510 1511 impl<T: SimdComplexField, D: Dim, S: RawStorageMut<T, D, D>> Matrix<T, D, D, S> { 1512 /// Sets `self` to its adjoint. 1513 #[deprecated(note = "Renamed to `self.adjoint_mut()`.")] conjugate_transform_mut(&mut self)1514 pub fn conjugate_transform_mut(&mut self) { 1515 self.adjoint_mut() 1516 } 1517 1518 /// Sets `self` to its adjoint (aka. conjugate-transpose). adjoint_mut(&mut self)1519 pub fn adjoint_mut(&mut self) { 1520 assert!( 1521 self.is_square(), 1522 "Unable to transpose a non-square matrix in-place." 1523 ); 1524 1525 let dim = self.shape().0; 1526 1527 for i in 0..dim { 1528 for j in 0..i { 1529 unsafe { 1530 let ref_ij = self.get_unchecked((i, j)).clone(); 1531 let ref_ji = self.get_unchecked((j, i)).clone(); 1532 let conj_ij = ref_ij.simd_conjugate(); 1533 let conj_ji = ref_ji.simd_conjugate(); 1534 *self.get_unchecked_mut((i, j)) = conj_ji; 1535 *self.get_unchecked_mut((j, i)) = conj_ij; 1536 } 1537 } 1538 1539 { 1540 let diag = unsafe { self.get_unchecked_mut((i, i)) }; 1541 *diag = diag.clone().simd_conjugate(); 1542 } 1543 } 1544 } 1545 } 1546 1547 impl<T: Scalar, D: Dim, S: RawStorage<T, D, D>> SquareMatrix<T, D, S> { 1548 /// The diagonal of this matrix. 1549 #[inline] 1550 #[must_use] diagonal(&self) -> OVector<T, D> where DefaultAllocator: Allocator<T, D>,1551 pub fn diagonal(&self) -> OVector<T, D> 1552 where 1553 DefaultAllocator: Allocator<T, D>, 1554 { 1555 self.map_diagonal(|e| e) 1556 } 1557 1558 /// Apply the given function to this matrix's diagonal and returns it. 1559 /// 1560 /// This is a more efficient version of `self.diagonal().map(f)` since this 1561 /// allocates only once. 1562 #[must_use] map_diagonal<T2: Scalar>(&self, mut f: impl FnMut(T) -> T2) -> OVector<T2, D> where DefaultAllocator: Allocator<T2, D>,1563 pub fn map_diagonal<T2: Scalar>(&self, mut f: impl FnMut(T) -> T2) -> OVector<T2, D> 1564 where 1565 DefaultAllocator: Allocator<T2, D>, 1566 { 1567 assert!( 1568 self.is_square(), 1569 "Unable to get the diagonal of a non-square matrix." 1570 ); 1571 1572 let dim = self.shape_generic().0; 1573 let mut res = Matrix::uninit(dim, Const::<1>); 1574 1575 for i in 0..dim.value() { 1576 // Safety: all indices are in range. 1577 unsafe { 1578 *res.vget_unchecked_mut(i) = 1579 MaybeUninit::new(f(self.get_unchecked((i, i)).clone())); 1580 } 1581 } 1582 1583 // Safety: res is now fully initialized. 1584 unsafe { res.assume_init() } 1585 } 1586 1587 /// Computes a trace of a square matrix, i.e., the sum of its diagonal elements. 1588 #[inline] 1589 #[must_use] trace(&self) -> T where T: Scalar + Zero + ClosedAdd,1590 pub fn trace(&self) -> T 1591 where 1592 T: Scalar + Zero + ClosedAdd, 1593 { 1594 assert!( 1595 self.is_square(), 1596 "Cannot compute the trace of non-square matrix." 1597 ); 1598 1599 let dim = self.shape_generic().0; 1600 let mut res = T::zero(); 1601 1602 for i in 0..dim.value() { 1603 res += unsafe { self.get_unchecked((i, i)).clone() }; 1604 } 1605 1606 res 1607 } 1608 } 1609 1610 impl<T: SimdComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> { 1611 /// The symmetric part of `self`, i.e., `0.5 * (self + self.transpose())`. 1612 #[inline] 1613 #[must_use] symmetric_part(&self) -> OMatrix<T, D, D> where DefaultAllocator: Allocator<T, D, D>,1614 pub fn symmetric_part(&self) -> OMatrix<T, D, D> 1615 where 1616 DefaultAllocator: Allocator<T, D, D>, 1617 { 1618 assert!( 1619 self.is_square(), 1620 "Cannot compute the symmetric part of a non-square matrix." 1621 ); 1622 let mut tr = self.transpose(); 1623 tr += self; 1624 tr *= crate::convert::<_, T>(0.5); 1625 tr 1626 } 1627 1628 /// The hermitian part of `self`, i.e., `0.5 * (self + self.adjoint())`. 1629 #[inline] 1630 #[must_use] hermitian_part(&self) -> OMatrix<T, D, D> where DefaultAllocator: Allocator<T, D, D>,1631 pub fn hermitian_part(&self) -> OMatrix<T, D, D> 1632 where 1633 DefaultAllocator: Allocator<T, D, D>, 1634 { 1635 assert!( 1636 self.is_square(), 1637 "Cannot compute the hermitian part of a non-square matrix." 1638 ); 1639 1640 let mut tr = self.adjoint(); 1641 tr += self; 1642 tr *= crate::convert::<_, T>(0.5); 1643 tr 1644 } 1645 } 1646 1647 impl<T: Scalar + Zero + One, D: DimAdd<U1> + IsNotStaticOne, S: RawStorage<T, D, D>> 1648 Matrix<T, D, D, S> 1649 { 1650 /// Yields the homogeneous matrix for this matrix, i.e., appending an additional dimension and 1651 /// and setting the diagonal element to `1`. 1652 #[inline] 1653 #[must_use] to_homogeneous(&self) -> OMatrix<T, DimSum<D, U1>, DimSum<D, U1>> where DefaultAllocator: Allocator<T, DimSum<D, U1>, DimSum<D, U1>>,1654 pub fn to_homogeneous(&self) -> OMatrix<T, DimSum<D, U1>, DimSum<D, U1>> 1655 where 1656 DefaultAllocator: Allocator<T, DimSum<D, U1>, DimSum<D, U1>>, 1657 { 1658 assert!( 1659 self.is_square(), 1660 "Only square matrices can currently be transformed to homogeneous coordinates." 1661 ); 1662 let dim = DimSum::<D, U1>::from_usize(self.nrows() + 1); 1663 let mut res = OMatrix::identity_generic(dim, dim); 1664 res.generic_slice_mut::<D, D>((0, 0), self.shape_generic()) 1665 .copy_from(self); 1666 res 1667 } 1668 } 1669 1670 impl<T: Scalar + Zero, D: DimAdd<U1>, S: RawStorage<T, D>> Vector<T, D, S> { 1671 /// Computes the coordinates in projective space of this vector, i.e., appends a `0` to its 1672 /// coordinates. 1673 #[inline] 1674 #[must_use] to_homogeneous(&self) -> OVector<T, DimSum<D, U1>> where DefaultAllocator: Allocator<T, DimSum<D, U1>>,1675 pub fn to_homogeneous(&self) -> OVector<T, DimSum<D, U1>> 1676 where 1677 DefaultAllocator: Allocator<T, DimSum<D, U1>>, 1678 { 1679 self.push(T::zero()) 1680 } 1681 1682 /// Constructs a vector from coordinates in projective space, i.e., removes a `0` at the end of 1683 /// `self`. Returns `None` if this last component is not zero. 1684 #[inline] from_homogeneous<SB>(v: Vector<T, DimSum<D, U1>, SB>) -> Option<OVector<T, D>> where SB: RawStorage<T, DimSum<D, U1>>, DefaultAllocator: Allocator<T, D>,1685 pub fn from_homogeneous<SB>(v: Vector<T, DimSum<D, U1>, SB>) -> Option<OVector<T, D>> 1686 where 1687 SB: RawStorage<T, DimSum<D, U1>>, 1688 DefaultAllocator: Allocator<T, D>, 1689 { 1690 if v[v.len() - 1].is_zero() { 1691 let nrows = D::from_usize(v.len() - 1); 1692 Some(v.generic_slice((0, 0), (nrows, Const::<1>)).into_owned()) 1693 } else { 1694 None 1695 } 1696 } 1697 } 1698 1699 impl<T: Scalar, D: DimAdd<U1>, S: RawStorage<T, D>> Vector<T, D, S> { 1700 /// Constructs a new vector of higher dimension by appending `element` to the end of `self`. 1701 #[inline] 1702 #[must_use] push(&self, element: T) -> OVector<T, DimSum<D, U1>> where DefaultAllocator: Allocator<T, DimSum<D, U1>>,1703 pub fn push(&self, element: T) -> OVector<T, DimSum<D, U1>> 1704 where 1705 DefaultAllocator: Allocator<T, DimSum<D, U1>>, 1706 { 1707 let len = self.len(); 1708 let hnrows = DimSum::<D, U1>::from_usize(len + 1); 1709 let mut res = Matrix::uninit(hnrows, Const::<1>); 1710 // This is basically a copy_from except that we warp the copied 1711 // values into MaybeUninit. 1712 res.generic_slice_mut((0, 0), self.shape_generic()) 1713 .zip_apply(self, |out, e| *out = MaybeUninit::new(e)); 1714 res[(len, 0)] = MaybeUninit::new(element); 1715 1716 // Safety: res has been fully initialized. 1717 unsafe { res.assume_init() } 1718 } 1719 } 1720 1721 impl<T, R: Dim, C: Dim, S> AbsDiffEq for Matrix<T, R, C, S> 1722 where 1723 T: Scalar + AbsDiffEq, 1724 S: RawStorage<T, R, C>, 1725 T::Epsilon: Clone, 1726 { 1727 type Epsilon = T::Epsilon; 1728 1729 #[inline] default_epsilon() -> Self::Epsilon1730 fn default_epsilon() -> Self::Epsilon { 1731 T::default_epsilon() 1732 } 1733 1734 #[inline] abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool1735 fn abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool { 1736 self.iter() 1737 .zip(other.iter()) 1738 .all(|(a, b)| a.abs_diff_eq(b, epsilon.clone())) 1739 } 1740 } 1741 1742 impl<T, R: Dim, C: Dim, S> RelativeEq for Matrix<T, R, C, S> 1743 where 1744 T: Scalar + RelativeEq, 1745 S: Storage<T, R, C>, 1746 T::Epsilon: Clone, 1747 { 1748 #[inline] default_max_relative() -> Self::Epsilon1749 fn default_max_relative() -> Self::Epsilon { 1750 T::default_max_relative() 1751 } 1752 1753 #[inline] relative_eq( &self, other: &Self, epsilon: Self::Epsilon, max_relative: Self::Epsilon, ) -> bool1754 fn relative_eq( 1755 &self, 1756 other: &Self, 1757 epsilon: Self::Epsilon, 1758 max_relative: Self::Epsilon, 1759 ) -> bool { 1760 self.relative_eq(other, epsilon, max_relative) 1761 } 1762 } 1763 1764 impl<T, R: Dim, C: Dim, S> UlpsEq for Matrix<T, R, C, S> 1765 where 1766 T: Scalar + UlpsEq, 1767 S: RawStorage<T, R, C>, 1768 T::Epsilon: Clone, 1769 { 1770 #[inline] default_max_ulps() -> u321771 fn default_max_ulps() -> u32 { 1772 T::default_max_ulps() 1773 } 1774 1775 #[inline] ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool1776 fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool { 1777 assert!(self.shape() == other.shape()); 1778 self.iter() 1779 .zip(other.iter()) 1780 .all(|(a, b)| a.ulps_eq(b, epsilon.clone(), max_ulps.clone())) 1781 } 1782 } 1783 1784 impl<T, R: Dim, C: Dim, S> PartialOrd for Matrix<T, R, C, S> 1785 where 1786 T: Scalar + PartialOrd, 1787 S: RawStorage<T, R, C>, 1788 { 1789 #[inline] partial_cmp(&self, other: &Self) -> Option<Ordering>1790 fn partial_cmp(&self, other: &Self) -> Option<Ordering> { 1791 if self.shape() != other.shape() { 1792 return None; 1793 } 1794 1795 if self.nrows() == 0 || self.ncols() == 0 { 1796 return Some(Ordering::Equal); 1797 } 1798 1799 let mut first_ord = unsafe { 1800 self.data 1801 .get_unchecked_linear(0) 1802 .partial_cmp(other.data.get_unchecked_linear(0)) 1803 }; 1804 1805 if let Some(first_ord) = first_ord.as_mut() { 1806 let mut it = self.iter().zip(other.iter()); 1807 let _ = it.next(); // Drop the first elements (we already tested it). 1808 1809 for (left, right) in it { 1810 if let Some(ord) = left.partial_cmp(right) { 1811 match ord { 1812 Ordering::Equal => { /* Does not change anything. */ } 1813 Ordering::Less => { 1814 if *first_ord == Ordering::Greater { 1815 return None; 1816 } 1817 *first_ord = ord 1818 } 1819 Ordering::Greater => { 1820 if *first_ord == Ordering::Less { 1821 return None; 1822 } 1823 *first_ord = ord 1824 } 1825 } 1826 } else { 1827 return None; 1828 } 1829 } 1830 } 1831 1832 first_ord 1833 } 1834 1835 #[inline] lt(&self, right: &Self) -> bool1836 fn lt(&self, right: &Self) -> bool { 1837 assert_eq!( 1838 self.shape(), 1839 right.shape(), 1840 "Matrix comparison error: dimensions mismatch." 1841 ); 1842 self.iter().zip(right.iter()).all(|(a, b)| a.lt(b)) 1843 } 1844 1845 #[inline] le(&self, right: &Self) -> bool1846 fn le(&self, right: &Self) -> bool { 1847 assert_eq!( 1848 self.shape(), 1849 right.shape(), 1850 "Matrix comparison error: dimensions mismatch." 1851 ); 1852 self.iter().zip(right.iter()).all(|(a, b)| a.le(b)) 1853 } 1854 1855 #[inline] gt(&self, right: &Self) -> bool1856 fn gt(&self, right: &Self) -> bool { 1857 assert_eq!( 1858 self.shape(), 1859 right.shape(), 1860 "Matrix comparison error: dimensions mismatch." 1861 ); 1862 self.iter().zip(right.iter()).all(|(a, b)| a.gt(b)) 1863 } 1864 1865 #[inline] ge(&self, right: &Self) -> bool1866 fn ge(&self, right: &Self) -> bool { 1867 assert_eq!( 1868 self.shape(), 1869 right.shape(), 1870 "Matrix comparison error: dimensions mismatch." 1871 ); 1872 self.iter().zip(right.iter()).all(|(a, b)| a.ge(b)) 1873 } 1874 } 1875 1876 impl<T, R: Dim, C: Dim, S> Eq for Matrix<T, R, C, S> 1877 where 1878 T: Scalar + Eq, 1879 S: RawStorage<T, R, C>, 1880 { 1881 } 1882 1883 impl<T, R, R2, C, C2, S, S2> PartialEq<Matrix<T, R2, C2, S2>> for Matrix<T, R, C, S> 1884 where 1885 T: Scalar + PartialEq, 1886 C: Dim, 1887 C2: Dim, 1888 R: Dim, 1889 R2: Dim, 1890 S: RawStorage<T, R, C>, 1891 S2: RawStorage<T, R2, C2>, 1892 { 1893 #[inline] eq(&self, right: &Matrix<T, R2, C2, S2>) -> bool1894 fn eq(&self, right: &Matrix<T, R2, C2, S2>) -> bool { 1895 self.shape() == right.shape() && self.iter().zip(right.iter()).all(|(l, r)| l == r) 1896 } 1897 } 1898 1899 macro_rules! impl_fmt { 1900 ($trait: path, $fmt_str_without_precision: expr, $fmt_str_with_precision: expr) => { 1901 impl<T, R: Dim, C: Dim, S> $trait for Matrix<T, R, C, S> 1902 where 1903 T: Scalar + $trait, 1904 S: RawStorage<T, R, C>, 1905 { 1906 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { 1907 #[cfg(feature = "std")] 1908 fn val_width<T: Scalar + $trait>(val: &T, f: &mut fmt::Formatter<'_>) -> usize { 1909 match f.precision() { 1910 Some(precision) => format!($fmt_str_with_precision, val, precision) 1911 .chars() 1912 .count(), 1913 None => format!($fmt_str_without_precision, val).chars().count(), 1914 } 1915 } 1916 1917 #[cfg(not(feature = "std"))] 1918 fn val_width<T: Scalar + $trait>(_: &T, _: &mut fmt::Formatter<'_>) -> usize { 1919 4 1920 } 1921 1922 let (nrows, ncols) = self.shape(); 1923 1924 if nrows == 0 || ncols == 0 { 1925 return write!(f, "[ ]"); 1926 } 1927 1928 let mut max_length = 0; 1929 1930 for i in 0..nrows { 1931 for j in 0..ncols { 1932 max_length = crate::max(max_length, val_width(&self[(i, j)], f)); 1933 } 1934 } 1935 1936 let max_length_with_space = max_length + 1; 1937 1938 writeln!(f)?; 1939 writeln!( 1940 f, 1941 " ┌ {:>width$} ┐", 1942 "", 1943 width = max_length_with_space * ncols - 1 1944 )?; 1945 1946 for i in 0..nrows { 1947 write!(f, " │")?; 1948 for j in 0..ncols { 1949 let number_length = val_width(&self[(i, j)], f) + 1; 1950 let pad = max_length_with_space - number_length; 1951 write!(f, " {:>thepad$}", "", thepad = pad)?; 1952 match f.precision() { 1953 Some(precision) => { 1954 write!(f, $fmt_str_with_precision, (*self)[(i, j)], precision)? 1955 } 1956 None => write!(f, $fmt_str_without_precision, (*self)[(i, j)])?, 1957 } 1958 } 1959 writeln!(f, " │")?; 1960 } 1961 1962 writeln!( 1963 f, 1964 " └ {:>width$} ┘", 1965 "", 1966 width = max_length_with_space * ncols - 1 1967 )?; 1968 writeln!(f) 1969 } 1970 } 1971 }; 1972 } 1973 impl_fmt!(fmt::Display, "{}", "{:.1$}"); 1974 impl_fmt!(fmt::LowerExp, "{:e}", "{:.1$e}"); 1975 impl_fmt!(fmt::UpperExp, "{:E}", "{:.1$E}"); 1976 impl_fmt!(fmt::Octal, "{:o}", "{:1$o}"); 1977 impl_fmt!(fmt::LowerHex, "{:x}", "{:1$x}"); 1978 impl_fmt!(fmt::UpperHex, "{:X}", "{:1$X}"); 1979 impl_fmt!(fmt::Binary, "{:b}", "{:.1$b}"); 1980 impl_fmt!(fmt::Pointer, "{:p}", "{:.1$p}"); 1981 1982 #[cfg(test)] 1983 mod tests { 1984 #[test] empty_display()1985 fn empty_display() { 1986 let vec: Vec<f64> = Vec::new(); 1987 let dvector = crate::DVector::from_vec(vec); 1988 assert_eq!(format!("{}", dvector), "[ ]") 1989 } 1990 1991 #[test] lower_exp()1992 fn lower_exp() { 1993 let test = crate::Matrix2::new(1e6, 2e5, 2e-5, 1.); 1994 assert_eq!( 1995 format!("{:e}", test), 1996 r" 1997 ┌ ┐ 1998 │ 1e6 2e5 │ 1999 │ 2e-5 1e0 │ 2000 └ ┘ 2001 2002 " 2003 ) 2004 } 2005 } 2006 2007 /// # Cross product 2008 impl<T: Scalar + ClosedAdd + ClosedSub + ClosedMul, R: Dim, C: Dim, S: RawStorage<T, R, C>> 2009 Matrix<T, R, C, S> 2010 { 2011 /// The perpendicular product between two 2D column vectors, i.e. `a.x * b.y - a.y * b.x`. 2012 #[inline] 2013 #[must_use] perp<R2, C2, SB>(&self, b: &Matrix<T, R2, C2, SB>) -> T where R2: Dim, C2: Dim, SB: RawStorage<T, R2, C2>, ShapeConstraint: SameNumberOfRows<R, U2> + SameNumberOfColumns<C, U1> + SameNumberOfRows<R2, U2> + SameNumberOfColumns<C2, U1>,2014 pub fn perp<R2, C2, SB>(&self, b: &Matrix<T, R2, C2, SB>) -> T 2015 where 2016 R2: Dim, 2017 C2: Dim, 2018 SB: RawStorage<T, R2, C2>, 2019 ShapeConstraint: SameNumberOfRows<R, U2> 2020 + SameNumberOfColumns<C, U1> 2021 + SameNumberOfRows<R2, U2> 2022 + SameNumberOfColumns<C2, U1>, 2023 { 2024 assert!( 2025 self.shape() == (2, 1), 2026 "2D perpendicular product requires (2, 1) vector but found {:?}", 2027 self.shape() 2028 ); 2029 2030 unsafe { 2031 self.get_unchecked((0, 0)).clone() * b.get_unchecked((1, 0)).clone() 2032 - self.get_unchecked((1, 0)).clone() * b.get_unchecked((0, 0)).clone() 2033 } 2034 } 2035 2036 // TODO: use specialization instead of an assertion. 2037 /// The 3D cross product between two vectors. 2038 /// 2039 /// Panics if the shape is not 3D vector. In the future, this will be implemented only for 2040 /// dynamically-sized matrices and statically-sized 3D matrices. 2041 #[inline] 2042 #[must_use] cross<R2, C2, SB>(&self, b: &Matrix<T, R2, C2, SB>) -> MatrixCross<T, R, C, R2, C2> where R2: Dim, C2: Dim, SB: RawStorage<T, R2, C2>, DefaultAllocator: SameShapeAllocator<T, R, C, R2, C2>, ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>,2043 pub fn cross<R2, C2, SB>(&self, b: &Matrix<T, R2, C2, SB>) -> MatrixCross<T, R, C, R2, C2> 2044 where 2045 R2: Dim, 2046 C2: Dim, 2047 SB: RawStorage<T, R2, C2>, 2048 DefaultAllocator: SameShapeAllocator<T, R, C, R2, C2>, 2049 ShapeConstraint: SameNumberOfRows<R, R2> + SameNumberOfColumns<C, C2>, 2050 { 2051 let shape = self.shape(); 2052 assert_eq!(shape, b.shape(), "Vector cross product dimension mismatch."); 2053 assert!( 2054 (shape.0 == 3 && shape.1 == 1) || (shape.0 == 1 && shape.1 == 3), 2055 "Vector cross product dimension mismatch: must be (3, 1) or (1, 3) but found {:?}.", 2056 shape 2057 ); 2058 2059 if shape.0 == 3 { 2060 unsafe { 2061 // TODO: soooo ugly! 2062 let nrows = SameShapeR::<R, R2>::from_usize(3); 2063 let ncols = SameShapeC::<C, C2>::from_usize(1); 2064 let mut res = Matrix::uninit(nrows, ncols); 2065 2066 let ax = self.get_unchecked((0, 0)); 2067 let ay = self.get_unchecked((1, 0)); 2068 let az = self.get_unchecked((2, 0)); 2069 2070 let bx = b.get_unchecked((0, 0)); 2071 let by = b.get_unchecked((1, 0)); 2072 let bz = b.get_unchecked((2, 0)); 2073 2074 *res.get_unchecked_mut((0, 0)) = 2075 MaybeUninit::new(ay.clone() * bz.clone() - az.clone() * by.clone()); 2076 *res.get_unchecked_mut((1, 0)) = 2077 MaybeUninit::new(az.clone() * bx.clone() - ax.clone() * bz.clone()); 2078 *res.get_unchecked_mut((2, 0)) = 2079 MaybeUninit::new(ax.clone() * by.clone() - ay.clone() * bx.clone()); 2080 2081 // Safety: res is now fully initialized. 2082 res.assume_init() 2083 } 2084 } else { 2085 unsafe { 2086 // TODO: ugly! 2087 let nrows = SameShapeR::<R, R2>::from_usize(1); 2088 let ncols = SameShapeC::<C, C2>::from_usize(3); 2089 let mut res = Matrix::uninit(nrows, ncols); 2090 2091 let ax = self.get_unchecked((0, 0)); 2092 let ay = self.get_unchecked((0, 1)); 2093 let az = self.get_unchecked((0, 2)); 2094 2095 let bx = b.get_unchecked((0, 0)); 2096 let by = b.get_unchecked((0, 1)); 2097 let bz = b.get_unchecked((0, 2)); 2098 2099 *res.get_unchecked_mut((0, 0)) = 2100 MaybeUninit::new(ay.clone() * bz.clone() - az.clone() * by.clone()); 2101 *res.get_unchecked_mut((0, 1)) = 2102 MaybeUninit::new(az.clone() * bx.clone() - ax.clone() * bz.clone()); 2103 *res.get_unchecked_mut((0, 2)) = 2104 MaybeUninit::new(ax.clone() * by.clone() - ay.clone() * bx.clone()); 2105 2106 // Safety: res is now fully initialized. 2107 res.assume_init() 2108 } 2109 } 2110 } 2111 } 2112 2113 impl<T: Scalar + Field, S: RawStorage<T, U3>> Vector<T, U3, S> { 2114 /// Computes the matrix `M` such that for all vector `v` we have `M * v == self.cross(&v)`. 2115 #[inline] 2116 #[must_use] cross_matrix(&self) -> OMatrix<T, U3, U3>2117 pub fn cross_matrix(&self) -> OMatrix<T, U3, U3> { 2118 OMatrix::<T, U3, U3>::new( 2119 T::zero(), 2120 -self[2].clone(), 2121 self[1].clone(), 2122 self[2].clone(), 2123 T::zero(), 2124 -self[0].clone(), 2125 -self[1].clone(), 2126 self[0].clone(), 2127 T::zero(), 2128 ) 2129 } 2130 } 2131 2132 impl<T: SimdComplexField, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> { 2133 /// The smallest angle between two vectors. 2134 #[inline] 2135 #[must_use] angle<R2: Dim, C2: Dim, SB>(&self, other: &Matrix<T, R2, C2, SB>) -> T::SimdRealField where SB: Storage<T, R2, C2>, ShapeConstraint: DimEq<R, R2> + DimEq<C, C2>,2136 pub fn angle<R2: Dim, C2: Dim, SB>(&self, other: &Matrix<T, R2, C2, SB>) -> T::SimdRealField 2137 where 2138 SB: Storage<T, R2, C2>, 2139 ShapeConstraint: DimEq<R, R2> + DimEq<C, C2>, 2140 { 2141 let prod = self.dotc(other); 2142 let n1 = self.norm(); 2143 let n2 = other.norm(); 2144 2145 if n1.is_zero() || n2.is_zero() { 2146 T::SimdRealField::zero() 2147 } else { 2148 let cang = prod.simd_real() / (n1 * n2); 2149 cang.simd_clamp(-T::SimdRealField::one(), T::SimdRealField::one()) 2150 .simd_acos() 2151 } 2152 } 2153 } 2154 2155 impl<T, R: Dim, C: Dim, S> AbsDiffEq for Unit<Matrix<T, R, C, S>> 2156 where 2157 T: Scalar + AbsDiffEq, 2158 S: RawStorage<T, R, C>, 2159 T::Epsilon: Clone, 2160 { 2161 type Epsilon = T::Epsilon; 2162 2163 #[inline] default_epsilon() -> Self::Epsilon2164 fn default_epsilon() -> Self::Epsilon { 2165 T::default_epsilon() 2166 } 2167 2168 #[inline] abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool2169 fn abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool { 2170 self.as_ref().abs_diff_eq(other.as_ref(), epsilon) 2171 } 2172 } 2173 2174 impl<T, R: Dim, C: Dim, S> RelativeEq for Unit<Matrix<T, R, C, S>> 2175 where 2176 T: Scalar + RelativeEq, 2177 S: Storage<T, R, C>, 2178 T::Epsilon: Clone, 2179 { 2180 #[inline] default_max_relative() -> Self::Epsilon2181 fn default_max_relative() -> Self::Epsilon { 2182 T::default_max_relative() 2183 } 2184 2185 #[inline] relative_eq( &self, other: &Self, epsilon: Self::Epsilon, max_relative: Self::Epsilon, ) -> bool2186 fn relative_eq( 2187 &self, 2188 other: &Self, 2189 epsilon: Self::Epsilon, 2190 max_relative: Self::Epsilon, 2191 ) -> bool { 2192 self.as_ref() 2193 .relative_eq(other.as_ref(), epsilon, max_relative) 2194 } 2195 } 2196 2197 impl<T, R: Dim, C: Dim, S> UlpsEq for Unit<Matrix<T, R, C, S>> 2198 where 2199 T: Scalar + UlpsEq, 2200 S: RawStorage<T, R, C>, 2201 T::Epsilon: Clone, 2202 { 2203 #[inline] default_max_ulps() -> u322204 fn default_max_ulps() -> u32 { 2205 T::default_max_ulps() 2206 } 2207 2208 #[inline] ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool2209 fn ulps_eq(&self, other: &Self, epsilon: Self::Epsilon, max_ulps: u32) -> bool { 2210 self.as_ref().ulps_eq(other.as_ref(), epsilon, max_ulps) 2211 } 2212 } 2213 2214 impl<T, R, C, S> Hash for Matrix<T, R, C, S> 2215 where 2216 T: Scalar + Hash, 2217 R: Dim, 2218 C: Dim, 2219 S: RawStorage<T, R, C>, 2220 { hash<H: Hasher>(&self, state: &mut H)2221 fn hash<H: Hasher>(&self, state: &mut H) { 2222 let (nrows, ncols) = self.shape(); 2223 (nrows, ncols).hash(state); 2224 2225 for j in 0..ncols { 2226 for i in 0..nrows { 2227 unsafe { 2228 self.get_unchecked((i, j)).hash(state); 2229 } 2230 } 2231 } 2232 } 2233 } 2234