diff --git a/nalgebra-lapack/src/cholesky.rs b/nalgebra-lapack/src/cholesky.rs
index fbbf6fd4..bc3515a5 100644
--- a/nalgebra-lapack/src/cholesky.rs
+++ b/nalgebra-lapack/src/cholesky.rs
@@ -80,6 +80,7 @@ where
     }
 
     /// Retrieves the lower-triangular factor of the cholesky decomposition.
+    #[must_use]
     pub fn l(&self) -> OMatrix<T, D, D> {
         let mut res = self.l.clone();
         res.fill_upper_triangle(Zero::zero(), 1);
@@ -91,6 +92,7 @@ where
     ///
     /// This is an allocation-less version of `self.l()`. The values of the strict upper-triangular
     /// part are garbage and should be ignored by further computations.
+    #[must_use]
     pub fn l_dirty(&self) -> &OMatrix<T, D, D> {
         &self.l
     }
diff --git a/nalgebra-lapack/src/eigen.rs b/nalgebra-lapack/src/eigen.rs
index cbe1555a..1bca79a5 100644
--- a/nalgebra-lapack/src/eigen.rs
+++ b/nalgebra-lapack/src/eigen.rs
@@ -302,6 +302,7 @@ where
 
     /// The determinant of the decomposed matrix.
     #[inline]
+    #[must_use]
     pub fn determinant(&self) -> T {
         let mut det = T::one();
         for e in self.eigenvalues.iter() {
diff --git a/nalgebra-lapack/src/hessenberg.rs b/nalgebra-lapack/src/hessenberg.rs
index b025a8bf..c5765022 100644
--- a/nalgebra-lapack/src/hessenberg.rs
+++ b/nalgebra-lapack/src/hessenberg.rs
@@ -89,6 +89,7 @@ where
 
     /// Computes the hessenberg matrix of this decomposition.
     #[inline]
+    #[must_use]
     pub fn h(&self) -> OMatrix<T, D, D> {
         let mut h = self.h.clone_owned();
         h.fill_lower_triangle(T::zero(), 2);
@@ -109,6 +110,7 @@ where
 
     /// Computes the unitary matrix `Q` of this decomposition.
     #[inline]
+    #[must_use]
     pub fn q(&self) -> OMatrix<T, D, D> {
         let n = self.h.nrows() as i32;
         let mut q = self.h.clone_owned();
diff --git a/nalgebra-lapack/src/lu.rs b/nalgebra-lapack/src/lu.rs
index 5ca5f143..2130fc7e 100644
--- a/nalgebra-lapack/src/lu.rs
+++ b/nalgebra-lapack/src/lu.rs
@@ -85,6 +85,7 @@ where
 
     /// Gets the lower-triangular matrix part of the decomposition.
     #[inline]
+    #[must_use]
     pub fn l(&self) -> OMatrix<T, R, DimMinimum<R, C>> {
         let (nrows, ncols) = self.lu.data.shape();
         let mut res = self.lu.columns_generic(0, nrows.min(ncols)).into_owned();
@@ -97,6 +98,7 @@ where
 
     /// Gets the upper-triangular matrix part of the decomposition.
     #[inline]
+    #[must_use]
     pub fn u(&self) -> OMatrix<T, DimMinimum<R, C>, C> {
         let (nrows, ncols) = self.lu.data.shape();
         let mut res = self.lu.rows_generic(0, nrows.min(ncols)).into_owned();
@@ -111,6 +113,7 @@ where
     /// Computing the permutation matrix explicitly is costly and usually not necessary.
     /// To permute rows of a matrix or vector, use the method `self.permute(...)` instead.
     #[inline]
+    #[must_use]
     pub fn p(&self) -> OMatrix<T, R, R> {
         let (dim, _) = self.lu.data.shape();
         let mut id = Matrix::identity_generic(dim, dim);
@@ -124,6 +127,7 @@ where
 
     /// Gets the LAPACK permutation indices.
     #[inline]
+    #[must_use]
     pub fn permutation_indices(&self) -> &OVector<i32, DimMinimum<R, C>> {
         &self.p
     }
diff --git a/nalgebra-lapack/src/qr.rs b/nalgebra-lapack/src/qr.rs
index b58c3f35..7b2d5df6 100644
--- a/nalgebra-lapack/src/qr.rs
+++ b/nalgebra-lapack/src/qr.rs
@@ -92,6 +92,7 @@ where
 
     /// Retrieves the upper trapezoidal submatrix `R` of this decomposition.
     #[inline]
+    #[must_use]
     pub fn r(&self) -> OMatrix<T, DimMinimum<R, C>, C> {
         let (nrows, ncols) = self.qr.data.shape();
         self.qr.rows_generic(0, nrows.min(ncols)).upper_triangle()
@@ -117,6 +118,7 @@ where
 
     /// Computes the orthogonal matrix `Q` of this decomposition.
     #[inline]
+    #[must_use]
     pub fn q(&self) -> OMatrix<T, R, DimMinimum<R, C>> {
         let (nrows, ncols) = self.qr.data.shape();
         let min_nrows_ncols = nrows.min(ncols);
diff --git a/nalgebra-lapack/src/schur.rs b/nalgebra-lapack/src/schur.rs
index a850f01d..3bee2635 100644
--- a/nalgebra-lapack/src/schur.rs
+++ b/nalgebra-lapack/src/schur.rs
@@ -138,6 +138,7 @@ where
     /// Computes the real eigenvalues of the decomposed matrix.
     ///
     /// Return `None` if some eigenvalues are complex.
+    #[must_use]
     pub fn eigenvalues(&self) -> Option<OVector<T, D>> {
         if self.im.iter().all(|e| e.is_zero()) {
             Some(self.re.clone())
@@ -147,6 +148,7 @@ where
     }
 
     /// Computes the complex eigenvalues of the decomposed matrix.
+    #[must_use]
     pub fn complex_eigenvalues(&self) -> OVector<Complex<T>, D>
     where
         DefaultAllocator: Allocator<Complex<T>, D>,
diff --git a/nalgebra-lapack/src/svd.rs b/nalgebra-lapack/src/svd.rs
index f6b907b6..3357e621 100644
--- a/nalgebra-lapack/src/svd.rs
+++ b/nalgebra-lapack/src/svd.rs
@@ -175,6 +175,7 @@ macro_rules! svd_impl(
             ///
             /// All singular value below epsilon will be set to zero instead of being inverted.
             #[inline]
+            #[must_use]
             pub fn pseudo_inverse(&self, epsilon: $t) -> OMatrix<$t, C, R> {
                 let nrows           = self.u.data.shape().0;
                 let ncols           = self.vt.data.shape().1;
@@ -207,6 +208,7 @@ macro_rules! svd_impl(
             /// This is the number of singular values that are not too small (i.e. greater than
             /// the given `epsilon`).
             #[inline]
+            #[must_use]
             pub fn rank(&self, epsilon: $t) -> usize {
                 let mut i = 0;
 
diff --git a/nalgebra-lapack/src/symmetric_eigen.rs b/nalgebra-lapack/src/symmetric_eigen.rs
index 0cd6bfc1..d276437e 100644
--- a/nalgebra-lapack/src/symmetric_eigen.rs
+++ b/nalgebra-lapack/src/symmetric_eigen.rs
@@ -138,6 +138,7 @@ where
 
     /// The determinant of the decomposed matrix.
     #[inline]
+    #[must_use]
     pub fn determinant(&self) -> T {
         let mut det = T::one();
         for e in self.eigenvalues.iter() {
@@ -150,6 +151,7 @@ where
     /// Rebuild the original matrix.
     ///
     /// This is useful if some of the eigenvalues have been manually modified.
+    #[must_use]
     pub fn recompose(&self) -> OMatrix<T, D, D> {
         let mut u_t = self.eigenvectors.clone();
         for i in 0..self.eigenvalues.len() {
diff --git a/nalgebra-sparse/src/coo.rs b/nalgebra-sparse/src/coo.rs
index 364083a9..caf74654 100644
--- a/nalgebra-sparse/src/coo.rs
+++ b/nalgebra-sparse/src/coo.rs
@@ -173,12 +173,14 @@ impl<T> CooMatrix<T> {
 
     /// The number of rows in the matrix.
     #[inline]
+    #[must_use]
     pub fn nrows(&self) -> usize {
         self.nrows
     }
 
     /// The number of columns in the matrix.
     #[inline]
+    #[must_use]
     pub fn ncols(&self) -> usize {
         self.ncols
     }
@@ -189,21 +191,25 @@ impl<T> CooMatrix<T> {
     /// entries, then it may have a different number of non-zeros as reported by `nnz()` compared
     /// to its CSR representation.
     #[inline]
+    #[must_use]
     pub fn nnz(&self) -> usize {
         self.values.len()
     }
 
     /// The row indices of the explicitly stored entries.
+    #[must_use]
     pub fn row_indices(&self) -> &[usize] {
         &self.row_indices
     }
 
     /// The column indices of the explicitly stored entries.
+    #[must_use]
     pub fn col_indices(&self) -> &[usize] {
         &self.col_indices
     }
 
     /// The values of the explicitly stored entries.
+    #[must_use]
     pub fn values(&self) -> &[T] {
         &self.values
     }
diff --git a/nalgebra-sparse/src/cs.rs b/nalgebra-sparse/src/cs.rs
index d6f9b229..cde0a3e2 100644
--- a/nalgebra-sparse/src/cs.rs
+++ b/nalgebra-sparse/src/cs.rs
@@ -32,11 +32,13 @@ impl<T> CsMatrix<T> {
     }
 
     #[inline]
+    #[must_use]
     pub fn pattern(&self) -> &SparsityPattern {
         &self.sparsity_pattern
     }
 
     #[inline]
+    #[must_use]
     pub fn values(&self) -> &[T] {
         &self.values
     }
@@ -48,6 +50,7 @@ impl<T> CsMatrix<T> {
 
     /// Returns the raw data represented as a tuple `(major_offsets, minor_indices, values)`.
     #[inline]
+    #[must_use]
     pub fn cs_data(&self) -> (&[usize], &[usize], &[T]) {
         let pattern = self.pattern();
         (
@@ -88,6 +91,7 @@ impl<T> CsMatrix<T> {
 
     /// Internal method for simplifying access to a lane's data
     #[inline]
+    #[must_use]
     pub fn get_index_range(&self, row_index: usize) -> Option<Range<usize>> {
         let row_begin = *self.sparsity_pattern.major_offsets().get(row_index)?;
         let row_end = *self.sparsity_pattern.major_offsets().get(row_index + 1)?;
@@ -111,6 +115,7 @@ impl<T> CsMatrix<T> {
 
     /// Returns an entry for the given major/minor indices, or `None` if the indices are out
     /// of bounds.
+    #[must_use]
     pub fn get_entry(&self, major_index: usize, minor_index: usize) -> Option<SparseEntry<T>> {
         let row_range = self.get_index_range(major_index)?;
         let (_, minor_indices, values) = self.cs_data();
@@ -139,6 +144,7 @@ impl<T> CsMatrix<T> {
         get_mut_entry_from_slices(minor_dim, minor_indices, values, minor_index)
     }
 
+    #[must_use]
     pub fn get_lane(&self, index: usize) -> Option<CsLane<T>> {
         let range = self.get_index_range(index)?;
         let (_, minor_indices, values) = self.cs_data();
@@ -150,6 +156,7 @@ impl<T> CsMatrix<T> {
     }
 
     #[inline]
+    #[must_use]
     pub fn get_lane_mut(&mut self, index: usize) -> Option<CsLaneMut<T>> {
         let range = self.get_index_range(index)?;
         let minor_dim = self.pattern().minor_dim();
@@ -172,6 +179,7 @@ impl<T> CsMatrix<T> {
     }
 
     #[inline]
+    #[must_use]
     pub fn filter<P>(&self, predicate: P) -> Self
     where
         T: Clone,
@@ -207,6 +215,7 @@ impl<T> CsMatrix<T> {
     }
 
     /// Returns the diagonal of the matrix as a sparse matrix.
+    #[must_use]
     pub fn diagonal_as_matrix(&self) -> Self
     where
         T: Clone,
@@ -372,26 +381,31 @@ macro_rules! impl_cs_lane_common_methods {
     ($name:ty) => {
         impl<'a, T> $name {
             #[inline]
+            #[must_use]
             pub fn minor_dim(&self) -> usize {
                 self.minor_dim
             }
 
             #[inline]
+            #[must_use]
             pub fn nnz(&self) -> usize {
                 self.minor_indices.len()
             }
 
             #[inline]
+            #[must_use]
             pub fn minor_indices(&self) -> &[usize] {
                 self.minor_indices
             }
 
             #[inline]
+            #[must_use]
             pub fn values(&self) -> &[T] {
                 self.values
             }
 
             #[inline]
+            #[must_use]
             pub fn get_entry(&self, global_col_index: usize) -> Option<SparseEntry<T>> {
                 get_entry_from_slices(
                     self.minor_dim,
@@ -416,6 +430,7 @@ impl<'a, T> CsLaneMut<'a, T> {
         (self.minor_indices, self.values)
     }
 
+    #[must_use]
     pub fn get_entry_mut(&mut self, global_minor_index: usize) -> Option<SparseEntryMut<T>> {
         get_mut_entry_from_slices(
             self.minor_dim,
diff --git a/nalgebra-sparse/src/csc.rs b/nalgebra-sparse/src/csc.rs
index 6a88c1d5..a691e1d4 100644
--- a/nalgebra-sparse/src/csc.rs
+++ b/nalgebra-sparse/src/csc.rs
@@ -192,12 +192,14 @@ impl<T> CscMatrix<T> {
 
     /// The number of rows in the matrix.
     #[inline]
+    #[must_use]
     pub fn nrows(&self) -> usize {
         self.cs.pattern().minor_dim()
     }
 
     /// The number of columns in the matrix.
     #[inline]
+    #[must_use]
     pub fn ncols(&self) -> usize {
         self.cs.pattern().major_dim()
     }
@@ -208,24 +210,28 @@ impl<T> CscMatrix<T> {
     /// number of algebraically zero entries in the matrix. Explicitly stored entries can still
     /// be zero. Corresponds to the number of entries in the sparsity pattern.
     #[inline]
+    #[must_use]
     pub fn nnz(&self) -> usize {
         self.pattern().nnz()
     }
 
     /// The column offsets defining part of the CSC format.
     #[inline]
+    #[must_use]
     pub fn col_offsets(&self) -> &[usize] {
         self.pattern().major_offsets()
     }
 
     /// The row indices defining part of the CSC format.
     #[inline]
+    #[must_use]
     pub fn row_indices(&self) -> &[usize] {
         self.pattern().minor_indices()
     }
 
     /// The non-zero values defining part of the CSC format.
     #[inline]
+    #[must_use]
     pub fn values(&self) -> &[T] {
         self.cs.values()
     }
@@ -298,6 +304,7 @@ impl<T> CscMatrix<T> {
     /// ------
     /// Panics if column index is out of bounds.
     #[inline]
+    #[must_use]
     pub fn col(&self, index: usize) -> CscCol<T> {
         self.get_col(index).expect("Row index must be in bounds")
     }
@@ -315,12 +322,14 @@ impl<T> CscMatrix<T> {
 
     /// Return the column at the given column index, or `None` if out of bounds.
     #[inline]
+    #[must_use]
     pub fn get_col(&self, index: usize) -> Option<CscCol<T>> {
         self.cs.get_lane(index).map(|lane| CscCol { lane })
     }
 
     /// Mutable column access for the given column index, or `None` if out of bounds.
     #[inline]
+    #[must_use]
     pub fn get_col_mut(&mut self, index: usize) -> Option<CscColMut<T>> {
         self.cs.get_lane_mut(index).map(|lane| CscColMut { lane })
     }
@@ -381,6 +390,7 @@ impl<T> CscMatrix<T> {
     }
 
     /// Returns a reference to the underlying sparsity pattern.
+    #[must_use]
     pub fn pattern(&self) -> &SparsityPattern {
         self.cs.pattern()
     }
@@ -397,6 +407,7 @@ impl<T> CscMatrix<T> {
     ///
     /// Each call to this function incurs the cost of a binary search among the explicitly
     /// stored row entries for the given column.
+    #[must_use]
     pub fn get_entry(&self, row_index: usize, col_index: usize) -> Option<SparseEntry<T>> {
         self.cs.get_entry(col_index, row_index)
     }
@@ -422,6 +433,7 @@ impl<T> CscMatrix<T> {
     /// Panics
     /// ------
     /// Panics if `row_index` or `col_index` is out of bounds.
+    #[must_use]
     pub fn index_entry(&self, row_index: usize, col_index: usize) -> SparseEntry<T> {
         self.get_entry(row_index, col_index)
             .expect("Out of bounds matrix indices encountered")
@@ -441,6 +453,7 @@ impl<T> CscMatrix<T> {
     }
 
     /// Returns a triplet of slices `(col_offsets, row_indices, values)` that make up the CSC data.
+    #[must_use]
     pub fn csc_data(&self) -> (&[usize], &[usize], &[T]) {
         self.cs.cs_data()
     }
@@ -453,6 +466,7 @@ impl<T> CscMatrix<T> {
 
     /// Creates a sparse matrix that contains only the explicit entries decided by the
     /// given predicate.
+    #[must_use]
     pub fn filter<P>(&self, predicate: P) -> Self
     where
         T: Clone,
@@ -470,6 +484,7 @@ impl<T> CscMatrix<T> {
     /// Returns a new matrix representing the upper triangular part of this matrix.
     ///
     /// The result includes the diagonal of the matrix.
+    #[must_use]
     pub fn upper_triangle(&self) -> Self
     where
         T: Clone,
@@ -480,6 +495,7 @@ impl<T> CscMatrix<T> {
     /// Returns a new matrix representing the lower triangular part of this matrix.
     ///
     /// The result includes the diagonal of the matrix.
+    #[must_use]
     pub fn lower_triangle(&self) -> Self
     where
         T: Clone,
@@ -488,6 +504,7 @@ impl<T> CscMatrix<T> {
     }
 
     /// Returns the diagonal of the matrix as a sparse matrix.
+    #[must_use]
     pub fn diagonal_as_csc(&self) -> Self
     where
         T: Clone,
@@ -498,6 +515,7 @@ impl<T> CscMatrix<T> {
     }
 
     /// Compute the transpose of the matrix.
+    #[must_use]
     pub fn transpose(&self) -> CscMatrix<T>
     where
         T: Scalar,
@@ -617,24 +635,28 @@ macro_rules! impl_csc_col_common_methods {
         impl<'a, T> $name {
             /// The number of global rows in the column.
             #[inline]
+            #[must_use]
             pub fn nrows(&self) -> usize {
                 self.lane.minor_dim()
             }
 
             /// The number of non-zeros in this column.
             #[inline]
+            #[must_use]
             pub fn nnz(&self) -> usize {
                 self.lane.nnz()
             }
 
             /// The row indices corresponding to explicitly stored entries in this column.
             #[inline]
+            #[must_use]
             pub fn row_indices(&self) -> &[usize] {
                 self.lane.minor_indices()
             }
 
             /// The values corresponding to explicitly stored entries in this column.
             #[inline]
+            #[must_use]
             pub fn values(&self) -> &[T] {
                 self.lane.values()
             }
@@ -643,6 +665,7 @@ macro_rules! impl_csc_col_common_methods {
             ///
             /// Each call to this function incurs the cost of a binary search among the explicitly
             /// stored row entries.
+            #[must_use]
             pub fn get_entry(&self, global_row_index: usize) -> Option<SparseEntry<T>> {
                 self.lane.get_entry(global_row_index)
             }
@@ -669,6 +692,7 @@ impl<'a, T> CscColMut<'a, T> {
     }
 
     /// Returns a mutable entry for the given global row index.
+    #[must_use]
     pub fn get_entry_mut(&mut self, global_row_index: usize) -> Option<SparseEntryMut<T>> {
         self.lane.get_entry_mut(global_row_index)
     }
diff --git a/nalgebra-sparse/src/csr.rs b/nalgebra-sparse/src/csr.rs
index ded189eb..c2f9b669 100644
--- a/nalgebra-sparse/src/csr.rs
+++ b/nalgebra-sparse/src/csr.rs
@@ -192,12 +192,14 @@ impl<T> CsrMatrix<T> {
 
     /// The number of rows in the matrix.
     #[inline]
+    #[must_use]
     pub fn nrows(&self) -> usize {
         self.cs.pattern().major_dim()
     }
 
     /// The number of columns in the matrix.
     #[inline]
+    #[must_use]
     pub fn ncols(&self) -> usize {
         self.cs.pattern().minor_dim()
     }
@@ -208,12 +210,14 @@ impl<T> CsrMatrix<T> {
     /// number of algebraically zero entries in the matrix. Explicitly stored entries can still
     /// be zero. Corresponds to the number of entries in the sparsity pattern.
     #[inline]
+    #[must_use]
     pub fn nnz(&self) -> usize {
         self.cs.pattern().nnz()
     }
 
     /// The row offsets defining part of the CSR format.
     #[inline]
+    #[must_use]
     pub fn row_offsets(&self) -> &[usize] {
         let (offsets, _, _) = self.cs.cs_data();
         offsets
@@ -221,6 +225,7 @@ impl<T> CsrMatrix<T> {
 
     /// The column indices defining part of the CSR format.
     #[inline]
+    #[must_use]
     pub fn col_indices(&self) -> &[usize] {
         let (_, indices, _) = self.cs.cs_data();
         indices
@@ -228,6 +233,7 @@ impl<T> CsrMatrix<T> {
 
     /// The non-zero values defining part of the CSR format.
     #[inline]
+    #[must_use]
     pub fn values(&self) -> &[T] {
         self.cs.values()
     }
@@ -300,6 +306,7 @@ impl<T> CsrMatrix<T> {
     /// ------
     /// Panics if row index is out of bounds.
     #[inline]
+    #[must_use]
     pub fn row(&self, index: usize) -> CsrRow<T> {
         self.get_row(index).expect("Row index must be in bounds")
     }
@@ -317,12 +324,14 @@ impl<T> CsrMatrix<T> {
 
     /// Return the row at the given row index, or `None` if out of bounds.
     #[inline]
+    #[must_use]
     pub fn get_row(&self, index: usize) -> Option<CsrRow<T>> {
         self.cs.get_lane(index).map(|lane| CsrRow { lane })
     }
 
     /// Mutable row access for the given row index, or `None` if out of bounds.
     #[inline]
+    #[must_use]
     pub fn get_row_mut(&mut self, index: usize) -> Option<CsrRowMut<T>> {
         self.cs.get_lane_mut(index).map(|lane| CsrRowMut { lane })
     }
@@ -383,6 +392,7 @@ impl<T> CsrMatrix<T> {
     }
 
     /// Returns a reference to the underlying sparsity pattern.
+    #[must_use]
     pub fn pattern(&self) -> &SparsityPattern {
         self.cs.pattern()
     }
@@ -399,6 +409,7 @@ impl<T> CsrMatrix<T> {
     ///
     /// Each call to this function incurs the cost of a binary search among the explicitly
     /// stored column entries for the given row.
+    #[must_use]
     pub fn get_entry(&self, row_index: usize, col_index: usize) -> Option<SparseEntry<T>> {
         self.cs.get_entry(row_index, col_index)
     }
@@ -424,6 +435,7 @@ impl<T> CsrMatrix<T> {
     /// Panics
     /// ------
     /// Panics if `row_index` or `col_index` is out of bounds.
+    #[must_use]
     pub fn index_entry(&self, row_index: usize, col_index: usize) -> SparseEntry<T> {
         self.get_entry(row_index, col_index)
             .expect("Out of bounds matrix indices encountered")
@@ -443,6 +455,7 @@ impl<T> CsrMatrix<T> {
     }
 
     /// Returns a triplet of slices `(row_offsets, col_indices, values)` that make up the CSR data.
+    #[must_use]
     pub fn csr_data(&self) -> (&[usize], &[usize], &[T]) {
         self.cs.cs_data()
     }
@@ -455,6 +468,7 @@ impl<T> CsrMatrix<T> {
 
     /// Creates a sparse matrix that contains only the explicit entries decided by the
     /// given predicate.
+    #[must_use]
     pub fn filter<P>(&self, predicate: P) -> Self
     where
         T: Clone,
@@ -470,6 +484,7 @@ impl<T> CsrMatrix<T> {
     /// Returns a new matrix representing the upper triangular part of this matrix.
     ///
     /// The result includes the diagonal of the matrix.
+    #[must_use]
     pub fn upper_triangle(&self) -> Self
     where
         T: Clone,
@@ -480,6 +495,7 @@ impl<T> CsrMatrix<T> {
     /// Returns a new matrix representing the lower triangular part of this matrix.
     ///
     /// The result includes the diagonal of the matrix.
+    #[must_use]
     pub fn lower_triangle(&self) -> Self
     where
         T: Clone,
@@ -488,6 +504,7 @@ impl<T> CsrMatrix<T> {
     }
 
     /// Returns the diagonal of the matrix as a sparse matrix.
+    #[must_use]
     pub fn diagonal_as_csr(&self) -> Self
     where
         T: Clone,
@@ -498,6 +515,7 @@ impl<T> CsrMatrix<T> {
     }
 
     /// Compute the transpose of the matrix.
+    #[must_use]
     pub fn transpose(&self) -> CsrMatrix<T>
     where
         T: Scalar,
@@ -617,24 +635,28 @@ macro_rules! impl_csr_row_common_methods {
         impl<'a, T> $name {
             /// The number of global columns in the row.
             #[inline]
+            #[must_use]
             pub fn ncols(&self) -> usize {
                 self.lane.minor_dim()
             }
 
             /// The number of non-zeros in this row.
             #[inline]
+            #[must_use]
             pub fn nnz(&self) -> usize {
                 self.lane.nnz()
             }
 
             /// The column indices corresponding to explicitly stored entries in this row.
             #[inline]
+            #[must_use]
             pub fn col_indices(&self) -> &[usize] {
                 self.lane.minor_indices()
             }
 
             /// The values corresponding to explicitly stored entries in this row.
             #[inline]
+            #[must_use]
             pub fn values(&self) -> &[T] {
                 self.lane.values()
             }
@@ -644,6 +666,7 @@ macro_rules! impl_csr_row_common_methods {
             /// Each call to this function incurs the cost of a binary search among the explicitly
             /// stored column entries.
             #[inline]
+            #[must_use]
             pub fn get_entry(&self, global_col_index: usize) -> Option<SparseEntry<T>> {
                 self.lane.get_entry(global_col_index)
             }
@@ -673,6 +696,7 @@ impl<'a, T> CsrRowMut<'a, T> {
 
     /// Returns a mutable entry for the given global column index.
     #[inline]
+    #[must_use]
     pub fn get_entry_mut(&mut self, global_col_index: usize) -> Option<SparseEntryMut<T>> {
         self.lane.get_entry_mut(global_col_index)
     }
diff --git a/nalgebra-sparse/src/factorization/cholesky.rs b/nalgebra-sparse/src/factorization/cholesky.rs
index a18761c9..0acc428d 100644
--- a/nalgebra-sparse/src/factorization/cholesky.rs
+++ b/nalgebra-sparse/src/factorization/cholesky.rs
@@ -42,6 +42,7 @@ impl CscSymbolicCholesky {
     }
 
     /// The pattern of the Cholesky factor `L`.
+    #[must_use]
     pub fn l_pattern(&self) -> &SparsityPattern {
         &self.l_pattern
     }
@@ -171,6 +172,7 @@ impl<T: RealField> CscCholesky<T> {
     }
 
     /// Returns a reference to the Cholesky factor `L`.
+    #[must_use]
     pub fn l(&self) -> &CscMatrix<T> {
         &self.l_factor
     }
@@ -260,6 +262,7 @@ impl<T: RealField> CscCholesky<T> {
     /// # Panics
     ///
     /// Panics if `B` is not square.
+    #[must_use = "Did you mean to use solve_mut()?"]
     pub fn solve<'a>(&'a self, b: impl Into<DMatrixSlice<'a, T>>) -> DMatrix<T> {
         let b = b.into();
         let mut output = b.clone_owned();
diff --git a/nalgebra-sparse/src/lib.rs b/nalgebra-sparse/src/lib.rs
index 4b96717c..6f28917d 100644
--- a/nalgebra-sparse/src/lib.rs
+++ b/nalgebra-sparse/src/lib.rs
@@ -170,6 +170,7 @@ pub struct SparseFormatError {
 
 impl SparseFormatError {
     /// The type of error.
+    #[must_use]
     pub fn kind(&self) -> &SparseFormatErrorKind {
         &self.kind
     }
diff --git a/nalgebra-sparse/src/ops/mod.rs b/nalgebra-sparse/src/ops/mod.rs
index a6a21fbc..50d25e63 100644
--- a/nalgebra-sparse/src/ops/mod.rs
+++ b/nalgebra-sparse/src/ops/mod.rs
@@ -140,11 +140,13 @@ pub enum Op<T> {
 
 impl<T> Op<T> {
     /// Returns a reference to the inner value that the operation applies to.
+    #[must_use]
     pub fn inner_ref(&self) -> &T {
         self.as_ref().into_inner()
     }
 
     /// Returns an `Op` applied to a reference of the inner value.
+    #[must_use]
     pub fn as_ref(&self) -> Op<&T> {
         match self {
             Op::NoOp(obj) => Op::NoOp(&obj),
diff --git a/nalgebra-sparse/src/ops/serial/mod.rs b/nalgebra-sparse/src/ops/serial/mod.rs
index 82285e86..4b0cc904 100644
--- a/nalgebra-sparse/src/ops/serial/mod.rs
+++ b/nalgebra-sparse/src/ops/serial/mod.rs
@@ -96,11 +96,13 @@ impl OperationError {
     }
 
     /// The operation error kind.
+    #[must_use]
     pub fn kind(&self) -> &OperationErrorKind {
         &self.error_kind
     }
 
     /// The underlying error message.
+    #[must_use]
     pub fn message(&self) -> &str {
         self.message.as_str()
     }
diff --git a/nalgebra-sparse/src/pattern.rs b/nalgebra-sparse/src/pattern.rs
index 08552d6a..2e490285 100644
--- a/nalgebra-sparse/src/pattern.rs
+++ b/nalgebra-sparse/src/pattern.rs
@@ -60,18 +60,21 @@ impl SparsityPattern {
 
     /// The offsets for the major dimension.
     #[inline]
+    #[must_use]
     pub fn major_offsets(&self) -> &[usize] {
         &self.major_offsets
     }
 
     /// The indices for the minor dimension.
     #[inline]
+    #[must_use]
     pub fn minor_indices(&self) -> &[usize] {
         &self.minor_indices
     }
 
     /// The number of major lanes in the pattern.
     #[inline]
+    #[must_use]
     pub fn major_dim(&self) -> usize {
         assert!(self.major_offsets.len() > 0);
         self.major_offsets.len() - 1
@@ -79,12 +82,14 @@ impl SparsityPattern {
 
     /// The number of minor lanes in the pattern.
     #[inline]
+    #[must_use]
     pub fn minor_dim(&self) -> usize {
         self.minor_dim
     }
 
     /// The number of "non-zeros", i.e. explicitly stored entries in the pattern.
     #[inline]
+    #[must_use]
     pub fn nnz(&self) -> usize {
         self.minor_indices.len()
     }
@@ -96,12 +101,14 @@ impl SparsityPattern {
     ///
     /// Panics if `major_index` is out of bounds.
     #[inline]
+    #[must_use]
     pub fn lane(&self, major_index: usize) -> &[usize] {
         self.get_lane(major_index).unwrap()
     }
 
     /// Get the lane at the given index, or `None` if out of bounds.
     #[inline]
+    #[must_use]
     pub fn get_lane(&self, major_index: usize) -> Option<&[usize]> {
         let offset_begin = *self.major_offsets().get(major_index)?;
         let offset_end = *self.major_offsets().get(major_index + 1)?;
@@ -197,6 +204,7 @@ impl SparsityPattern {
     /// assert_eq!(entries, vec![(0, 0), (0, 2), (1, 1), (2, 0)]);
     /// ```
     ///
+    #[must_use]
     pub fn entries(&self) -> SparsityPatternIter {
         SparsityPatternIter::from_pattern(self)
     }
@@ -228,6 +236,7 @@ impl SparsityPattern {
     ///
     /// This is analogous to matrix transposition, i.e. an entry `(i, j)` becomes `(j, i)` in the
     /// new pattern.
+    #[must_use]
     pub fn transpose(&self) -> Self {
         // By using unit () values, we can use the same routines as for CSR/CSC matrices
         let values = vec![(); self.nnz()];
diff --git a/src/base/blas.rs b/src/base/blas.rs
index 05705f17..9617e46a 100644
--- a/src/base/blas.rs
+++ b/src/base/blas.rs
@@ -193,6 +193,7 @@ where
     /// ```
     ///
     #[inline]
+    #[must_use]
     pub fn dot<R2: Dim, C2: Dim, SB>(&self, rhs: &Matrix<T, R2, C2, SB>) -> T
     where
         SB: Storage<T, R2, C2>,
@@ -221,6 +222,7 @@ where
     /// assert_ne!(vec1.dotc(&vec2), vec1.dot(&vec2));
     /// ```
     #[inline]
+    #[must_use]
     pub fn dotc<R2: Dim, C2: Dim, SB>(&self, rhs: &Matrix<T, R2, C2, SB>) -> T
     where
         T: SimdComplexField,
@@ -248,6 +250,7 @@ where
     /// assert_eq!(mat1.tr_dot(&mat2), 9.1);
     /// ```
     #[inline]
+    #[must_use]
     pub fn tr_dot<R2: Dim, C2: Dim, SB>(&self, rhs: &Matrix<T, R2, C2, SB>) -> T
     where
         SB: Storage<T, R2, C2>,
diff --git a/src/base/componentwise.rs b/src/base/componentwise.rs
index 2ba49b2a..02b2cae6 100644
--- a/src/base/componentwise.rs
+++ b/src/base/componentwise.rs
@@ -28,6 +28,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(a.abs(), Matrix2::new(0.0, 1.0, 2.0, 3.0))
     /// ```
     #[inline]
+    #[must_use]
     pub fn abs(&self) -> OMatrix<T, R, C>
     where
         T: Signed,
@@ -49,6 +50,7 @@ macro_rules! component_binop_impl(
     ($($binop: ident, $binop_mut: ident, $binop_assign: ident, $cmpy: ident, $Trait: ident . $op: ident . $op_assign: ident, $desc:expr, $desc_cmpy:expr, $desc_mut:expr);* $(;)*) => {$(
         #[doc = $desc]
         #[inline]
+        #[must_use]
         pub fn $binop<R2, C2, SB>(&self, rhs: &Matrix<T, R2, C2, SB>) -> MatrixComponentOp<T, R1, C1, R2, C2>
             where T: $Trait,
                   R2: Dim, C2: Dim,
@@ -251,6 +253,7 @@ impl<T: Scalar, R1: Dim, C1: Dim, SA: Storage<T, R1, C1>> Matrix<T, R1, C1, SA>
     /// assert_eq!(u.inf(&v), expected)
     /// ```
     #[inline]
+    #[must_use]
     pub fn inf(&self, other: &Self) -> OMatrix<T, R1, C1>
     where
         T: SimdPartialOrd,
@@ -271,6 +274,7 @@ impl<T: Scalar, R1: Dim, C1: Dim, SA: Storage<T, R1, C1>> Matrix<T, R1, C1, SA>
     /// assert_eq!(u.sup(&v), expected)
     /// ```
     #[inline]
+    #[must_use]
     pub fn sup(&self, other: &Self) -> OMatrix<T, R1, C1>
     where
         T: SimdPartialOrd,
@@ -291,6 +295,7 @@ impl<T: Scalar, R1: Dim, C1: Dim, SA: Storage<T, R1, C1>> Matrix<T, R1, C1, SA>
     /// assert_eq!(u.inf_sup(&v), expected)
     /// ```
     #[inline]
+    #[must_use]
     pub fn inf_sup(&self, other: &Self) -> (OMatrix<T, R1, C1>, OMatrix<T, R1, C1>)
     where
         T: SimdPartialOrd,
diff --git a/src/base/edition.rs b/src/base/edition.rs
index 17c08c98..82c537fc 100644
--- a/src/base/edition.rs
+++ b/src/base/edition.rs
@@ -18,6 +18,7 @@ use crate::base::{DefaultAllocator, Matrix, OMatrix, RowVector, Scalar, Vector};
 impl<T: Scalar + Zero, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// Extracts the upper triangular part of this matrix (including the diagonal).
     #[inline]
+    #[must_use]
     pub fn upper_triangle(&self) -> OMatrix<T, R, C>
     where
         DefaultAllocator: Allocator<T, R, C>,
@@ -30,6 +31,7 @@ impl<T: Scalar + Zero, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 
     /// Extracts the lower triangular part of this matrix (including the diagonal).
     #[inline]
+    #[must_use]
     pub fn lower_triangle(&self) -> OMatrix<T, R, C>
     where
         DefaultAllocator: Allocator<T, R, C>,
@@ -42,6 +44,7 @@ impl<T: Scalar + Zero, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 
     /// Creates a new matrix by extracting the given set of rows from `self`.
     #[cfg(any(feature = "std", feature = "alloc"))]
+    #[must_use]
     pub fn select_rows<'a, I>(&self, irows: I) -> OMatrix<T, Dynamic, C>
     where
         I: IntoIterator<Item = &'a usize>,
@@ -78,6 +81,7 @@ impl<T: Scalar + Zero, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 
     /// Creates a new matrix by extracting the given set of columns from `self`.
     #[cfg(any(feature = "std", feature = "alloc"))]
+    #[must_use]
     pub fn select_columns<'a, I>(&self, icols: I) -> OMatrix<T, R, Dynamic>
     where
         I: IntoIterator<Item = &'a usize>,
diff --git a/src/base/indexing.rs b/src/base/indexing.rs
index 3e0c085e..5bc06020 100644
--- a/src/base/indexing.rs
+++ b/src/base/indexing.rs
@@ -485,6 +485,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// Produces a view of the data at the given index, or
     /// `None` if the index is out of bounds.
     #[inline]
+    #[must_use]
     pub fn get<'a, I>(&'a self, index: I) -> Option<I::Output>
     where
         I: MatrixIndex<'a, T, R, C, S>,
@@ -495,6 +496,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// Produces a mutable view of the data at the given index, or
     /// `None` if the index is out of bounds.
     #[inline]
+    #[must_use]
     pub fn get_mut<'a, I>(&'a mut self, index: I) -> Option<I::OutputMut>
     where
         S: StorageMut<T, R, C>,
@@ -506,6 +508,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// Produces a view of the data at the given index, or
     /// panics if the index is out of bounds.
     #[inline]
+    #[must_use]
     pub fn index<'a, I>(&'a self, index: I) -> I::Output
     where
         I: MatrixIndex<'a, T, R, C, S>,
@@ -527,6 +530,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// Produces a view of the data at the given index, without doing
     /// any bounds checking.
     #[inline]
+    #[must_use]
     pub unsafe fn get_unchecked<'a, I>(&'a self, index: I) -> I::Output
     where
         I: MatrixIndex<'a, T, R, C, S>,
@@ -537,6 +541,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// Returns a mutable view of the data at the given index, without doing
     /// any bounds checking.
     #[inline]
+    #[must_use]
     pub unsafe fn get_unchecked_mut<'a, I>(&'a mut self, index: I) -> I::OutputMut
     where
         S: StorageMut<T, R, C>,
diff --git a/src/base/interpolation.rs b/src/base/interpolation.rs
index c008743a..d5661e40 100644
--- a/src/base/interpolation.rs
+++ b/src/base/interpolation.rs
@@ -20,6 +20,7 @@ impl<T: Scalar + Zero + One + ClosedAdd + ClosedSub + ClosedMul, D: Dim, S: Stor
     /// let y = Vector3::new(10.0, 20.0, 30.0);
     /// assert_eq!(x.lerp(&y, 0.1), Vector3::new(1.9, 3.8, 5.7));
     /// ```
+    #[must_use]
     pub fn lerp<S2: Storage<T, D>>(&self, rhs: &Vector<T, D, S2>, t: T) -> OVector<T, D>
     where
         DefaultAllocator: Allocator<T, D>,
@@ -45,6 +46,7 @@ impl<T: Scalar + Zero + One + ClosedAdd + ClosedSub + ClosedMul, D: Dim, S: Stor
     ///
     /// assert_eq!(v, v2.normalize());
     /// ```
+    #[must_use]
     pub fn slerp<S2: Storage<T, D>>(&self, rhs: &Vector<T, D, S2>, t: T) -> OVector<T, D>
     where
         T: RealField,
@@ -72,6 +74,7 @@ impl<T: RealField, D: Dim, S: Storage<T, D>> Unit<Vector<T, D, S>> {
     ///
     /// assert_eq!(v, v2);
     /// ```
+    #[must_use]
     pub fn slerp<S2: Storage<T, D>>(
         &self,
         rhs: &Unit<Vector<T, D, S2>>,
@@ -89,6 +92,7 @@ impl<T: RealField, D: Dim, S: Storage<T, D>> Unit<Vector<T, D, S>> {
     ///
     /// Returns `None` if the two vectors are almost collinear and with opposite direction
     /// (in this case, there is an infinity of possible results).
+    #[must_use]
     pub fn try_slerp<S2: Storage<T, D>>(
         &self,
         rhs: &Unit<Vector<T, D, S2>>,
diff --git a/src/base/matrix.rs b/src/base/matrix.rs
index bc9942be..5214bb6f 100644
--- a/src/base/matrix.rs
+++ b/src/base/matrix.rs
@@ -441,6 +441,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// let mat = Matrix3x4::<f32>::zeros();
     /// assert_eq!(mat.shape(), (3, 4));
     #[inline]
+    #[must_use]
     pub fn shape(&self) -> (usize, usize) {
         let (nrows, ncols) = self.data.shape();
         (nrows.value(), ncols.value())
@@ -455,6 +456,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// let mat = Matrix3x4::<f32>::zeros();
     /// assert_eq!(mat.nrows(), 3);
     #[inline]
+    #[must_use]
     pub fn nrows(&self) -> usize {
         self.shape().0
     }
@@ -468,6 +470,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// let mat = Matrix3x4::<f32>::zeros();
     /// assert_eq!(mat.ncols(), 4);
     #[inline]
+    #[must_use]
     pub fn ncols(&self) -> usize {
         self.shape().1
     }
@@ -483,6 +486,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// // The column strides is the number of steps (here 2) multiplied by the corresponding dimension.
     /// assert_eq!(mat.strides(), (1, 10));
     #[inline]
+    #[must_use]
     pub fn strides(&self) -> (usize, usize) {
         let (srows, scols) = self.data.strides();
         (srows.value(), scols.value())
@@ -501,6 +505,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(m[i], m[3]);
     /// ```
     #[inline]
+    #[must_use]
     pub fn vector_to_matrix_index(&self, i: usize) -> (usize, usize) {
         let (nrows, ncols) = self.shape();
 
@@ -529,6 +534,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(unsafe { *ptr }, m[0]);
     /// ```
     #[inline]
+    #[must_use]
     pub fn as_ptr(&self) -> *const T {
         self.data.ptr()
     }
@@ -537,6 +543,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     ///
     /// See `relative_eq` from the `RelativeEq` trait for more details.
     #[inline]
+    #[must_use]
     pub fn relative_eq<R2, C2, SB>(
         &self,
         other: &Matrix<T, R2, C2, SB>,
@@ -559,6 +566,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 
     /// Tests whether `self` and `rhs` are exactly equal.
     #[inline]
+    #[must_use]
     pub fn eq<R2, C2, SB>(&self, other: &Matrix<T, R2, C2, SB>) -> bool
     where
         T: PartialEq,
@@ -609,6 +617,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 
     /// Clones this matrix to one that owns its data.
     #[inline]
+    #[must_use]
     pub fn clone_owned(&self) -> OMatrix<T, R, C>
     where
         DefaultAllocator: Allocator<T, R, C>,
@@ -619,6 +628,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// Clones this matrix into one that owns its data. The actual type of the result depends on
     /// matrix storage combination rules for addition.
     #[inline]
+    #[must_use]
     pub fn clone_owned_sum<R2, C2>(&self) -> MatrixSum<T, R, C, R2, C2>
     where
         R2: Dim,
@@ -692,6 +702,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// Returns a matrix containing the result of `f` applied to each of its entries.
     #[inline]
+    #[must_use]
     pub fn map<T2: Scalar, F: FnMut(T) -> T2>(&self, mut f: F) -> OMatrix<T2, R, C>
     where
         DefaultAllocator: Allocator<T2, R, C>,
@@ -738,6 +749,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// - If the matrix has has least one component, then `init_f` is called with the first component
     /// to compute the initial value. Folding then continues on all the remaining components of the matrix.
     #[inline]
+    #[must_use]
     pub fn fold_with<T2>(
         &self,
         init_f: impl FnOnce(Option<&T>) -> T2,
@@ -751,6 +763,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// Returns a matrix containing the result of `f` applied to each of its entries. Unlike `map`,
     /// `f` also gets passed the row and column index, i.e. `f(row, col, value)`.
     #[inline]
+    #[must_use]
     pub fn map_with_location<T2: Scalar, F: FnMut(usize, usize, T) -> T2>(
         &self,
         mut f: F,
@@ -778,6 +791,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// Returns a matrix containing the result of `f` applied to each entries of `self` and
     /// `rhs`.
     #[inline]
+    #[must_use]
     pub fn zip_map<T2, N3, S2, F>(&self, rhs: &Matrix<T2, R, C, S2>, mut f: F) -> OMatrix<N3, R, C>
     where
         T2: Scalar,
@@ -813,6 +827,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// Returns a matrix containing the result of `f` applied to each entries of `self` and
     /// `b`, and `c`.
     #[inline]
+    #[must_use]
     pub fn zip_zip_map<T2, N3, N4, S2, S3, F>(
         &self,
         b: &Matrix<T2, R, C, S2>,
@@ -860,6 +875,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 
     /// Folds a function `f` on each entry of `self`.
     #[inline]
+    #[must_use]
     pub fn fold<Acc>(&self, init: Acc, mut f: impl FnMut(Acc, T) -> Acc) -> Acc {
         let (nrows, ncols) = self.data.shape();
 
@@ -879,6 +895,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 
     /// Folds a function `f` on each pairs of entries from `self` and `rhs`.
     #[inline]
+    #[must_use]
     pub fn zip_fold<T2, R2, C2, S2, Acc>(
         &self,
         rhs: &Matrix<T2, R2, C2, S2>,
@@ -1238,6 +1255,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: StorageMut<T, R, C>> Matrix<T, R, C, S> {
 impl<T: Scalar, D: Dim, S: Storage<T, D>> Vector<T, D, S> {
     /// Gets a reference to the i-th element of this column vector without bound checking.
     #[inline]
+    #[must_use]
     pub unsafe fn vget_unchecked(&self, i: usize) -> &T {
         debug_assert!(i < self.nrows(), "Vector index out of bounds.");
         let i = i * self.strides().0;
@@ -1248,6 +1266,7 @@ impl<T: Scalar, D: Dim, S: Storage<T, D>> Vector<T, D, S> {
 impl<T: Scalar, D: Dim, S: StorageMut<T, D>> Vector<T, D, S> {
     /// Gets a mutable reference to the i-th element of this column vector without bound checking.
     #[inline]
+    #[must_use]
     pub unsafe fn vget_unchecked_mut(&mut self, i: usize) -> &mut T {
         debug_assert!(i < self.nrows(), "Vector index out of bounds.");
         let i = i * self.strides().0;
@@ -1258,6 +1277,7 @@ impl<T: Scalar, D: Dim, S: StorageMut<T, D>> Vector<T, D, S> {
 impl<T: Scalar, R: Dim, C: Dim, S: ContiguousStorage<T, R, C>> Matrix<T, R, C, S> {
     /// Extracts a slice containing the entire matrix entries ordered column-by-columns.
     #[inline]
+    #[must_use]
     pub fn as_slice(&self) -> &[T] {
         self.data.as_slice()
     }
@@ -1266,6 +1286,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: ContiguousStorage<T, R, C>> Matrix<T, R, C, S
 impl<T: Scalar, R: Dim, C: Dim, S: ContiguousStorageMut<T, R, C>> Matrix<T, R, C, S> {
     /// Extracts a mutable slice containing the entire matrix entries ordered column-by-columns.
     #[inline]
+    #[must_use]
     pub fn as_mut_slice(&mut self) -> &mut [T] {
         self.data.as_mut_slice()
     }
@@ -1446,6 +1467,7 @@ impl<T: SimdComplexField, D: Dim, S: StorageMut<T, D, D>> Matrix<T, D, D, S> {
 impl<T: Scalar, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
     /// The diagonal of this matrix.
     #[inline]
+    #[must_use]
     pub fn diagonal(&self) -> OVector<T, D>
     where
         DefaultAllocator: Allocator<T, D>,
@@ -1457,6 +1479,7 @@ impl<T: Scalar, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
     ///
     /// This is a more efficient version of `self.diagonal().map(f)` since this
     /// allocates only once.
+    #[must_use]
     pub fn map_diagonal<T2: Scalar>(&self, mut f: impl FnMut(T) -> T2) -> OVector<T2, D>
     where
         DefaultAllocator: Allocator<T2, D>,
@@ -1481,6 +1504,7 @@ impl<T: Scalar, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
 
     /// Computes a trace of a square matrix, i.e., the sum of its diagonal elements.
     #[inline]
+    #[must_use]
     pub fn trace(&self) -> T
     where
         T: Scalar + Zero + ClosedAdd,
@@ -1504,6 +1528,7 @@ impl<T: Scalar, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
 impl<T: SimdComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
     /// The symmetric part of `self`, i.e., `0.5 * (self + self.transpose())`.
     #[inline]
+    #[must_use]
     pub fn symmetric_part(&self) -> OMatrix<T, D, D>
     where
         DefaultAllocator: Allocator<T, D, D>,
@@ -1520,6 +1545,7 @@ impl<T: SimdComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
 
     /// The hermitian part of `self`, i.e., `0.5 * (self + self.adjoint())`.
     #[inline]
+    #[must_use]
     pub fn hermitian_part(&self) -> OMatrix<T, D, D>
     where
         DefaultAllocator: Allocator<T, D, D>,
@@ -1542,6 +1568,7 @@ impl<T: Scalar + Zero + One, D: DimAdd<U1> + IsNotStaticOne, S: Storage<T, D, D>
     /// Yields the homogeneous matrix for this matrix, i.e., appending an additional dimension and
     /// and setting the diagonal element to `1`.
     #[inline]
+    #[must_use]
     pub fn to_homogeneous(&self) -> OMatrix<T, DimSum<D, U1>, DimSum<D, U1>>
     where
         DefaultAllocator: Allocator<T, DimSum<D, U1>, DimSum<D, U1>>,
@@ -1562,6 +1589,7 @@ impl<T: Scalar + Zero, D: DimAdd<U1>, S: Storage<T, D>> Vector<T, D, S> {
     /// Computes the coordinates in projective space of this vector, i.e., appends a `0` to its
     /// coordinates.
     #[inline]
+    #[must_use]
     pub fn to_homogeneous(&self) -> OVector<T, DimSum<D, U1>>
     where
         DefaultAllocator: Allocator<T, DimSum<D, U1>>,
@@ -1589,6 +1617,7 @@ impl<T: Scalar + Zero, D: DimAdd<U1>, S: Storage<T, D>> Vector<T, D, S> {
 impl<T: Scalar + Zero, D: DimAdd<U1>, S: Storage<T, D>> Vector<T, D, S> {
     /// Constructs a new vector of higher dimension by appending `element` to the end of `self`.
     #[inline]
+    #[must_use]
     pub fn push(&self, element: T) -> OVector<T, DimSum<D, U1>>
     where
         DefaultAllocator: Allocator<T, DimSum<D, U1>>,
@@ -1891,6 +1920,7 @@ impl<T: Scalar + ClosedAdd + ClosedSub + ClosedMul, R: Dim, C: Dim, S: Storage<T
 {
     /// The perpendicular product between two 2D column vectors, i.e. `a.x * b.y - a.y * b.x`.
     #[inline]
+    #[must_use]
     pub fn perp<R2, C2, SB>(&self, b: &Matrix<T, R2, C2, SB>) -> T
     where
         R2: Dim,
@@ -1920,6 +1950,7 @@ impl<T: Scalar + ClosedAdd + ClosedSub + ClosedMul, R: Dim, C: Dim, S: Storage<T
     /// Panics if the shape is not 3D vector. In the future, this will be implemented only for
     /// dynamically-sized matrices and statically-sized 3D matrices.
     #[inline]
+    #[must_use]
     pub fn cross<R2, C2, SB>(&self, b: &Matrix<T, R2, C2, SB>) -> MatrixCross<T, R, C, R2, C2>
     where
         R2: Dim,
@@ -1993,6 +2024,7 @@ impl<T: Scalar + ClosedAdd + ClosedSub + ClosedMul, R: Dim, C: Dim, S: Storage<T
 impl<T: Scalar + Field, S: Storage<T, U3>> Vector<T, U3, S> {
     /// Computes the matrix `M` such that for all vector `v` we have `M * v == self.cross(&v)`.
     #[inline]
+    #[must_use]
     pub fn cross_matrix(&self) -> OMatrix<T, U3, U3> {
         OMatrix::<T, U3, U3>::new(
             T::zero(),
@@ -2011,6 +2043,7 @@ impl<T: Scalar + Field, S: Storage<T, U3>> Vector<T, U3, S> {
 impl<T: SimdComplexField, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// The smallest angle between two vectors.
     #[inline]
+    #[must_use]
     pub fn angle<R2: Dim, C2: Dim, SB>(&self, other: &Matrix<T, R2, C2, SB>) -> T::SimdRealField
     where
         SB: Storage<T, R2, C2>,
diff --git a/src/base/matrix_slice.rs b/src/base/matrix_slice.rs
index 6b03c7b2..ffe95c8a 100644
--- a/src/base/matrix_slice.rs
+++ b/src/base/matrix_slice.rs
@@ -812,6 +812,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// Slices a sub-matrix containing the rows indexed by the range `rows` and the columns indexed
     /// by the range `cols`.
     #[inline]
+    #[must_use]
     pub fn slice_range<RowRange, ColRange>(
         &self,
         rows: RowRange,
@@ -830,6 +831,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 
     /// Slice containing all the rows indexed by the range `rows`.
     #[inline]
+    #[must_use]
     pub fn rows_range<RowRange: SliceRange<R>>(
         &self,
         rows: RowRange,
@@ -839,6 +841,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 
     /// Slice containing all the columns indexed by the range `rows`.
     #[inline]
+    #[must_use]
     pub fn columns_range<ColRange: SliceRange<C>>(
         &self,
         cols: ColRange,
diff --git a/src/base/min_max.rs b/src/base/min_max.rs
index d532d87a..83e62d10 100644
--- a/src/base/min_max.rs
+++ b/src/base/min_max.rs
@@ -13,6 +13,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(Vector3::new(-1.0, -2.0, -3.0).amax(), 3.0);
     /// ```
     #[inline]
+    #[must_use]
     pub fn amax(&self) -> T
     where
         T: Zero + SimdSigned + SimdPartialOrd,
@@ -33,6 +34,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     ///     Complex::new(1.0, 3.0)).camax(), 5.0);
     /// ```
     #[inline]
+    #[must_use]
     pub fn camax(&self) -> T::SimdRealField
     where
         T: SimdComplexField,
@@ -52,6 +54,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(Vector3::new(5u32, 2, 3).max(), 5);
     /// ```
     #[inline]
+    #[must_use]
     pub fn max(&self) -> T
     where
         T: SimdPartialOrd + Zero,
@@ -70,6 +73,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(Vector3::new(10.0, 2.0, 30.0).amin(), 2.0);
     /// ```
     #[inline]
+    #[must_use]
     pub fn amin(&self) -> T
     where
         T: Zero + SimdPartialOrd + SimdSigned,
@@ -90,6 +94,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     ///     Complex::new(1.0, 3.0)).camin(), 3.0);
     /// ```
     #[inline]
+    #[must_use]
     pub fn camin(&self) -> T::SimdRealField
     where
         T: SimdComplexField,
@@ -112,6 +117,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(Vector3::new(5u32, 2, 3).min(), 2);
     /// ```
     #[inline]
+    #[must_use]
     pub fn min(&self) -> T
     where
         T: SimdPartialOrd + Zero,
@@ -136,6 +142,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(mat.icamax_full(), (1, 0));
     /// ```
     #[inline]
+    #[must_use]
     pub fn icamax_full(&self) -> (usize, usize)
     where
         T: ComplexField,
@@ -172,6 +179,7 @@ impl<T: Scalar + PartialOrd + Signed, R: Dim, C: Dim, S: Storage<T, R, C>> Matri
     /// assert_eq!(mat.iamax_full(), (1, 2));
     /// ```
     #[inline]
+    #[must_use]
     pub fn iamax_full(&self) -> (usize, usize) {
         assert!(!self.is_empty(), "The input matrix must not be empty.");
 
@@ -209,6 +217,7 @@ impl<T: Scalar, D: Dim, S: Storage<T, D>> Vector<T, D, S> {
     /// assert_eq!(vec.icamax(), 2);
     /// ```
     #[inline]
+    #[must_use]
     pub fn icamax(&self) -> usize
     where
         T: ComplexField,
@@ -240,6 +249,7 @@ impl<T: Scalar, D: Dim, S: Storage<T, D>> Vector<T, D, S> {
     /// assert_eq!(vec.argmax(), (2, 13));
     /// ```
     #[inline]
+    #[must_use]
     pub fn argmax(&self) -> (usize, T)
     where
         T: PartialOrd,
@@ -271,6 +281,7 @@ impl<T: Scalar, D: Dim, S: Storage<T, D>> Vector<T, D, S> {
     /// assert_eq!(vec.imax(), 2);
     /// ```
     #[inline]
+    #[must_use]
     pub fn imax(&self) -> usize
     where
         T: PartialOrd,
@@ -288,6 +299,7 @@ impl<T: Scalar, D: Dim, S: Storage<T, D>> Vector<T, D, S> {
     /// assert_eq!(vec.iamax(), 1);
     /// ```
     #[inline]
+    #[must_use]
     pub fn iamax(&self) -> usize
     where
         T: PartialOrd + Signed,
@@ -319,6 +331,7 @@ impl<T: Scalar, D: Dim, S: Storage<T, D>> Vector<T, D, S> {
     /// assert_eq!(vec.argmin(), (1, -15));
     /// ```
     #[inline]
+    #[must_use]
     pub fn argmin(&self) -> (usize, T)
     where
         T: PartialOrd,
@@ -350,6 +363,7 @@ impl<T: Scalar, D: Dim, S: Storage<T, D>> Vector<T, D, S> {
     /// assert_eq!(vec.imin(), 1);
     /// ```
     #[inline]
+    #[must_use]
     pub fn imin(&self) -> usize
     where
         T: PartialOrd,
@@ -367,6 +381,7 @@ impl<T: Scalar, D: Dim, S: Storage<T, D>> Vector<T, D, S> {
     /// assert_eq!(vec.iamin(), 0);
     /// ```
     #[inline]
+    #[must_use]
     pub fn iamin(&self) -> usize
     where
         T: PartialOrd + Signed,
diff --git a/src/base/norm.rs b/src/base/norm.rs
index 251bf661..09e11f7e 100644
--- a/src/base/norm.rs
+++ b/src/base/norm.rs
@@ -158,6 +158,7 @@ impl<T: SimdComplexField> Norm<T> for UniformNorm {
 impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// The squared L2 norm of this vector.
     #[inline]
+    #[must_use]
     pub fn norm_squared(&self) -> T::SimdRealField
     where
         T: SimdComplexField,
@@ -176,6 +177,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     ///
     /// Use `.apply_norm` to apply a custom norm.
     #[inline]
+    #[must_use]
     pub fn norm(&self) -> T::SimdRealField
     where
         T: SimdComplexField,
@@ -187,6 +189,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     ///
     /// Use `.apply_metric_distance` to apply a custom norm.
     #[inline]
+    #[must_use]
     pub fn metric_distance<R2, C2, S2>(&self, rhs: &Matrix<T, R2, C2, S2>) -> T::SimdRealField
     where
         T: SimdComplexField,
@@ -211,6 +214,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(v.apply_norm(&EuclideanNorm), v.norm());
     /// ```
     #[inline]
+    #[must_use]
     pub fn apply_norm(&self, norm: &impl Norm<T>) -> T::SimdRealField
     where
         T: SimdComplexField,
@@ -233,6 +237,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(v1.apply_metric_distance(&v2, &EuclideanNorm), (v1 - v2).norm());
     /// ```
     #[inline]
+    #[must_use]
     pub fn apply_metric_distance<R2, C2, S2>(
         &self,
         rhs: &Matrix<T, R2, C2, S2>,
@@ -254,6 +259,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     ///
     /// This function is simply implemented as a call to `norm()`
     #[inline]
+    #[must_use]
     pub fn magnitude(&self) -> T::SimdRealField
     where
         T: SimdComplexField,
@@ -267,6 +273,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     ///
     /// This function is simply implemented as a call to `norm_squared()`
     #[inline]
+    #[must_use]
     pub fn magnitude_squared(&self) -> T::SimdRealField
     where
         T: SimdComplexField,
@@ -298,6 +305,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 
     /// The Lp norm of this matrix.
     #[inline]
+    #[must_use]
     pub fn lp_norm(&self, p: i32) -> T::SimdRealField
     where
         T: SimdComplexField,
@@ -341,6 +349,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 
     /// Returns a new vector with the same magnitude as `self` clamped between `0.0` and `max`.
     #[inline]
+    #[must_use]
     pub fn cap_magnitude(&self, max: T::RealField) -> OMatrix<T, R, C>
     where
         T: ComplexField,
@@ -357,6 +366,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 
     /// Returns a new vector with the same magnitude as `self` clamped between `0.0` and `max`.
     #[inline]
+    #[must_use]
     pub fn simd_cap_magnitude(&self, max: T::SimdRealField) -> OMatrix<T, R, C>
     where
         T: SimdComplexField,
diff --git a/src/base/ops.rs b/src/base/ops.rs
index 8917efb1..d9be71a7 100644
--- a/src/base/ops.rs
+++ b/src/base/ops.rs
@@ -676,6 +676,7 @@ where
 {
     /// Equivalent to `self.transpose() * rhs`.
     #[inline]
+    #[must_use]
     pub fn tr_mul<R2: Dim, C2: Dim, SB>(&self, rhs: &Matrix<T, R2, C2, SB>) -> OMatrix<T, C1, C2>
     where
         SB: Storage<T, R2, C2>,
@@ -692,6 +693,7 @@ where
 
     /// Equivalent to `self.adjoint() * rhs`.
     #[inline]
+    #[must_use]
     pub fn ad_mul<R2: Dim, C2: Dim, SB>(&self, rhs: &Matrix<T, R2, C2, SB>) -> OMatrix<T, C1, C2>
     where
         T: SimdComplexField,
@@ -801,6 +803,7 @@ where
 
     /// The kronecker product of two matrices (aka. tensor product of the corresponding linear
     /// maps).
+    #[must_use]
     pub fn kronecker<R2: Dim, C2: Dim, SB>(
         &self,
         rhs: &Matrix<T, R2, C2, SB>,
diff --git a/src/base/properties.rs b/src/base/properties.rs
index 1cbff719..9e250119 100644
--- a/src/base/properties.rs
+++ b/src/base/properties.rs
@@ -20,6 +20,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(mat.len(), 12);
     /// ```
     #[inline]
+    #[must_use]
     pub fn len(&self) -> usize {
         let (nrows, ncols) = self.shape();
         nrows * ncols
@@ -35,12 +36,14 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert!(!mat.is_empty());
     /// ```
     #[inline]
+    #[must_use]
     pub fn is_empty(&self) -> bool {
         self.len() == 0
     }
 
     /// Indicates if this is a square matrix.
     #[inline]
+    #[must_use]
     pub fn is_square(&self) -> bool {
         let (nrows, ncols) = self.shape();
         nrows == ncols
@@ -52,6 +55,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// If the matrix is diagonal, this checks that diagonal elements (i.e. at coordinates `(i, i)`
     /// for i from `0` to `min(R, C)`) are equal one; and that all other elements are zero.
     #[inline]
+    #[must_use]
     pub fn is_identity(&self, eps: T::Epsilon) -> bool
     where
         T: Zero + One + RelativeEq,
@@ -112,6 +116,7 @@ impl<T: ComplexField, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// In this definition `Id` is approximately equal to the identity matrix with a relative error
     /// equal to `eps`.
     #[inline]
+    #[must_use]
     pub fn is_orthogonal(&self, eps: T::Epsilon) -> bool
     where
         T: Zero + One + ClosedAdd + ClosedMul + RelativeEq,
@@ -129,6 +134,7 @@ where
 {
     /// Checks that this matrix is orthogonal and has a determinant equal to 1.
     #[inline]
+    #[must_use]
     pub fn is_special_orthogonal(&self, eps: T) -> bool
     where
         D: DimMin<D, Output = D>,
@@ -139,6 +145,7 @@ where
 
     /// Returns `true` if this matrix is invertible.
     #[inline]
+    #[must_use]
     pub fn is_invertible(&self) -> bool {
         // TODO: improve this?
         self.clone_owned().try_inverse().is_some()
diff --git a/src/base/statistics.rs b/src/base/statistics.rs
index 8f6c77ce..59d78482 100644
--- a/src/base/statistics.rs
+++ b/src/base/statistics.rs
@@ -9,6 +9,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// Returns a row vector where each element is the result of the application of `f` on the
     /// corresponding column of the original matrix.
     #[inline]
+    #[must_use]
     pub fn compress_rows(
         &self,
         f: impl Fn(VectorSlice<T, R, S::RStride, S::CStride>) -> T,
@@ -35,6 +36,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     ///
     /// This is the same as `self.compress_rows(f).transpose()`.
     #[inline]
+    #[must_use]
     pub fn compress_rows_tr(
         &self,
         f: impl Fn(VectorSlice<T, R, S::RStride, S::CStride>) -> T,
@@ -58,6 +60,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
 
     /// Returns a column vector resulting from the folding of `f` on each column of this matrix.
     #[inline]
+    #[must_use]
     pub fn compress_columns(
         &self,
         init: OVector<T, R>,
@@ -95,6 +98,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(m.sum(), 21.0);
     /// ```
     #[inline]
+    #[must_use]
     pub fn sum(&self) -> T
     where
         T: ClosedAdd + Zero,
@@ -120,6 +124,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(mint.row_sum(), RowVector2::new(9,12));
     /// ```
     #[inline]
+    #[must_use]
     pub fn row_sum(&self) -> RowOVector<T, C>
     where
         T: ClosedAdd + Zero,
@@ -144,6 +149,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(mint.row_sum_tr(), Vector2::new(9,12));
     /// ```
     #[inline]
+    #[must_use]
     pub fn row_sum_tr(&self) -> OVector<T, C>
     where
         T: ClosedAdd + Zero,
@@ -168,6 +174,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(mint.column_sum(), Vector3::new(3,7,11));
     /// ```
     #[inline]
+    #[must_use]
     pub fn column_sum(&self) -> OVector<T, R>
     where
         T: ClosedAdd + Zero,
@@ -197,6 +204,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_relative_eq!(m.variance(), 35.0 / 12.0, epsilon = 1.0e-8);
     /// ```
     #[inline]
+    #[must_use]
     pub fn variance(&self) -> T
     where
         T: Field + SupersetOf<f64>,
@@ -226,6 +234,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(m.row_variance(), RowVector3::new(2.25, 2.25, 2.25));
     /// ```
     #[inline]
+    #[must_use]
     pub fn row_variance(&self) -> RowOVector<T, C>
     where
         T: Field + SupersetOf<f64>,
@@ -246,6 +255,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(m.row_variance_tr(), Vector3::new(2.25, 2.25, 2.25));
     /// ```
     #[inline]
+    #[must_use]
     pub fn row_variance_tr(&self) -> OVector<T, C>
     where
         T: Field + SupersetOf<f64>,
@@ -267,6 +277,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_relative_eq!(m.column_variance(), Vector2::new(2.0 / 3.0, 2.0 / 3.0), epsilon = 1.0e-8);
     /// ```
     #[inline]
+    #[must_use]
     pub fn column_variance(&self) -> OVector<T, R>
     where
         T: Field + SupersetOf<f64>,
@@ -306,6 +317,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(m.mean(), 3.5);
     /// ```
     #[inline]
+    #[must_use]
     pub fn mean(&self) -> T
     where
         T: Field + SupersetOf<f64>,
@@ -331,6 +343,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(m.row_mean(), RowVector3::new(2.5, 3.5, 4.5));
     /// ```
     #[inline]
+    #[must_use]
     pub fn row_mean(&self) -> RowOVector<T, C>
     where
         T: Field + SupersetOf<f64>,
@@ -351,6 +364,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(m.row_mean_tr(), Vector3::new(2.5, 3.5, 4.5));
     /// ```
     #[inline]
+    #[must_use]
     pub fn row_mean_tr(&self) -> OVector<T, C>
     where
         T: Field + SupersetOf<f64>,
@@ -371,6 +385,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: Storage<T, R, C>> Matrix<T, R, C, S> {
     /// assert_eq!(m.column_mean(), Vector2::new(2.0, 5.0));
     /// ```
     #[inline]
+    #[must_use]
     pub fn column_mean(&self) -> OVector<T, R>
     where
         T: Field + SupersetOf<f64>,
diff --git a/src/base/swizzle.rs b/src/base/swizzle.rs
index 1e2f0e4c..25d6375f 100644
--- a/src/base/swizzle.rs
+++ b/src/base/swizzle.rs
@@ -8,6 +8,7 @@ macro_rules! impl_swizzle {
             $(
                 /// Builds a new vector from components of `self`.
                 #[inline]
+                #[must_use]
                 pub fn $name(&self) -> $Result<T>
                 where D::Typenum: Cmp<typenum::$BaseDim, Output=Greater> {
                     $Result::new($(self[$i].inlined_clone()),*)
diff --git a/src/base/vec_storage.rs b/src/base/vec_storage.rs
index 8b2d2ef2..14490674 100644
--- a/src/base/vec_storage.rs
+++ b/src/base/vec_storage.rs
@@ -95,6 +95,7 @@ impl<T, R: Dim, C: Dim> VecStorage<T, R, C> {
 
     /// The underlying data storage.
     #[inline]
+    #[must_use]
     pub fn as_vec(&self) -> &Vec<T> {
         &self.data
     }
@@ -129,12 +130,14 @@ impl<T, R: Dim, C: Dim> VecStorage<T, R, C> {
 
     /// The number of elements on the underlying vector.
     #[inline]
+    #[must_use]
     pub fn len(&self) -> usize {
         self.data.len()
     }
 
     /// Returns true if the underlying vector contains no elements.
     #[inline]
+    #[must_use]
     pub fn is_empty(&self) -> bool {
         self.len() == 0
     }
diff --git a/src/geometry/dual_quaternion.rs b/src/geometry/dual_quaternion.rs
index b6b127ce..83912b9e 100644
--- a/src/geometry/dual_quaternion.rs
+++ b/src/geometry/dual_quaternion.rs
@@ -232,6 +232,7 @@ where
     /// ));
     /// ```
     #[inline]
+    #[must_use]
     pub fn lerp(&self, other: &Self, t: T) -> Self {
         self * (T::one() - t) + other * t
     }
@@ -381,6 +382,7 @@ where
     /// ));
     /// ```
     #[inline]
+    #[must_use]
     pub fn dual_quaternion(&self) -> &DualQuaternion<T> {
         self.as_ref()
     }
@@ -463,7 +465,6 @@ where
     /// assert_relative_eq!(inv * unit, UnitDualQuaternion::identity(), epsilon = 1.0e-6);
     /// ```
     #[inline]
-    #[must_use = "Did you mean to use inverse_mut()?"]
     pub fn inverse_mut(&mut self) {
         let quat = self.as_mut_unchecked();
         quat.real = Unit::new_unchecked(quat.real).inverse().into_inner();
@@ -486,6 +487,7 @@ where
     /// assert_relative_eq!(dq_to * dq1, dq2, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn isometry_to(&self, other: &Self) -> Self {
         other / self
     }
@@ -518,6 +520,7 @@ where
     /// );
     /// ```
     #[inline]
+    #[must_use]
     pub fn lerp(&self, other: &Self, t: T) -> DualQuaternion<T> {
         self.as_ref().lerp(other.as_ref(), t)
     }
@@ -546,6 +549,7 @@ where
     /// ), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn nlerp(&self, other: &Self, t: T) -> Self {
         let mut res = self.lerp(other, t);
         let _ = res.normalize_mut();
@@ -581,6 +585,7 @@ where
     /// );
     /// assert_relative_eq!(dq.translation().vector.y, 3.0, epsilon = 1.0e-6);
     #[inline]
+    #[must_use]
     pub fn sclerp(&self, other: &Self, t: T) -> Self
     where
         T: RealField,
@@ -600,6 +605,7 @@ where
     /// * `epsilon`: the value below which the sinus of the angle separating both quaternion
     /// must be to return `None`.
     #[inline]
+    #[must_use]
     pub fn try_sclerp(&self, other: &Self, t: T, epsilon: T) -> Option<Self>
     where
         T: RealField,
@@ -667,6 +673,7 @@ where
     /// );
     /// ```
     #[inline]
+    #[must_use]
     pub fn rotation(&self) -> UnitQuaternion<T> {
         Unit::new_unchecked(self.as_ref().real)
     }
@@ -686,6 +693,7 @@ where
     /// );
     /// ```
     #[inline]
+    #[must_use]
     pub fn translation(&self) -> Translation3<T> {
         let two = T::one() + T::one();
         Translation3::from(
@@ -712,6 +720,7 @@ where
     /// assert_relative_eq!(iso.translation.vector, translation, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn to_isometry(&self) -> Isometry3<T> {
         Isometry3::from_parts(self.translation(), self.rotation())
     }
@@ -735,6 +744,7 @@ where
     /// );
     /// ```
     #[inline]
+    #[must_use]
     pub fn transform_point(&self, pt: &Point3<T>) -> Point3<T> {
         self * pt
     }
@@ -758,6 +768,7 @@ where
     /// );
     /// ```
     #[inline]
+    #[must_use]
     pub fn transform_vector(&self, v: &Vector3<T>) -> Vector3<T> {
         self * v
     }
@@ -781,6 +792,7 @@ where
     /// );
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_point(&self, pt: &Point3<T>) -> Point3<T> {
         self.inverse() * pt
     }
@@ -805,6 +817,7 @@ where
     /// );
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_vector(&self, v: &Vector3<T>) -> Vector3<T> {
         self.inverse() * v
     }
@@ -830,6 +843,7 @@ where
     /// );
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_unit_vector(&self, v: &Unit<Vector3<T>>) -> Unit<Vector3<T>> {
         self.inverse() * v
     }
@@ -857,6 +871,7 @@ where
     /// assert_relative_eq!(dq.to_homogeneous(), expected, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn to_homogeneous(&self) -> Matrix4<T> {
         self.to_isometry().to_homogeneous()
     }
diff --git a/src/geometry/isometry.rs b/src/geometry/isometry.rs
index d307c393..6d3e9127 100755
--- a/src/geometry/isometry.rs
+++ b/src/geometry/isometry.rs
@@ -267,6 +267,7 @@ where
     /// assert_eq!(iso1.inverse() * iso2, iso1.inv_mul(&iso2));
     /// ```
     #[inline]
+    #[must_use]
     pub fn inv_mul(&self, rhs: &Isometry<T, R, D>) -> Self {
         let inv_rot1 = self.rotation.inverse();
         let tr_12 = rhs.translation.vector.clone() - self.translation.vector.clone();
@@ -384,6 +385,7 @@ where
     /// assert_relative_eq!(transformed_point, Point3::new(3.0, 2.0, 2.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn transform_point(&self, pt: &Point<T, D>) -> Point<T, D> {
         self * pt
     }
@@ -407,6 +409,7 @@ where
     /// assert_relative_eq!(transformed_point, Vector3::new(3.0, 2.0, -1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn transform_vector(&self, v: &SVector<T, D>) -> SVector<T, D> {
         self * v
     }
@@ -429,6 +432,7 @@ where
     /// assert_relative_eq!(transformed_point, Point3::new(0.0, 2.0, 1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_point(&self, pt: &Point<T, D>) -> Point<T, D> {
         self.rotation
             .inverse_transform_point(&(pt - &self.translation.vector))
@@ -453,6 +457,7 @@ where
     /// assert_relative_eq!(transformed_point, Vector3::new(-3.0, 2.0, 1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_vector(&self, v: &SVector<T, D>) -> SVector<T, D> {
         self.rotation.inverse_transform_vector(v)
     }
@@ -476,6 +481,7 @@ where
     /// assert_relative_eq!(transformed_point, -Vector3::y_axis(), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_unit_vector(&self, v: &Unit<SVector<T, D>>) -> Unit<SVector<T, D>> {
         self.rotation.inverse_transform_unit_vector(v)
     }
@@ -505,6 +511,7 @@ impl<T: SimdRealField, R, const D: usize> Isometry<T, R, D> {
     /// assert_relative_eq!(iso.to_homogeneous(), expected, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn to_homogeneous(&self) -> OMatrix<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>
     where
         Const<D>: DimNameAdd<U1>,
@@ -536,6 +543,7 @@ impl<T: SimdRealField, R, const D: usize> Isometry<T, R, D> {
     /// assert_relative_eq!(iso.to_matrix(), expected, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn to_matrix(&self) -> OMatrix<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>
     where
         Const<D>: DimNameAdd<U1>,
diff --git a/src/geometry/isometry_interpolation.rs b/src/geometry/isometry_interpolation.rs
index 2ee5461a..356dbdad 100644
--- a/src/geometry/isometry_interpolation.rs
+++ b/src/geometry/isometry_interpolation.rs
@@ -26,6 +26,7 @@ impl<T: SimdRealField> Isometry3<T> {
     /// assert_eq!(iso3.rotation.euler_angles(), (std::f32::consts::FRAC_PI_2, 0.0, 0.0));
     /// ```
     #[inline]
+    #[must_use]
     pub fn lerp_slerp(&self, other: &Self, t: T) -> Self
     where
         T: RealField,
@@ -59,6 +60,7 @@ impl<T: SimdRealField> Isometry3<T> {
     /// assert_eq!(iso3.rotation.euler_angles(), (std::f32::consts::FRAC_PI_2, 0.0, 0.0));
     /// ```
     #[inline]
+    #[must_use]
     pub fn try_lerp_slerp(&self, other: &Self, t: T, epsilon: T) -> Option<Self>
     where
         T: RealField,
@@ -94,6 +96,7 @@ impl<T: SimdRealField> IsometryMatrix3<T> {
     /// assert_eq!(iso3.rotation.euler_angles(), (std::f32::consts::FRAC_PI_2, 0.0, 0.0));
     /// ```
     #[inline]
+    #[must_use]
     pub fn lerp_slerp(&self, other: &Self, t: T) -> Self
     where
         T: RealField,
@@ -127,6 +130,7 @@ impl<T: SimdRealField> IsometryMatrix3<T> {
     /// assert_eq!(iso3.rotation.euler_angles(), (std::f32::consts::FRAC_PI_2, 0.0, 0.0));
     /// ```
     #[inline]
+    #[must_use]
     pub fn try_lerp_slerp(&self, other: &Self, t: T, epsilon: T) -> Option<Self>
     where
         T: RealField,
@@ -163,6 +167,7 @@ impl<T: SimdRealField> Isometry2<T> {
     /// assert_relative_eq!(iso3.rotation.angle(), std::f32::consts::FRAC_PI_2);
     /// ```
     #[inline]
+    #[must_use]
     pub fn lerp_slerp(&self, other: &Self, t: T) -> Self
     where
         T: RealField,
@@ -199,6 +204,7 @@ impl<T: SimdRealField> IsometryMatrix2<T> {
     /// assert_relative_eq!(iso3.rotation.angle(), std::f32::consts::FRAC_PI_2);
     /// ```
     #[inline]
+    #[must_use]
     pub fn lerp_slerp(&self, other: &Self, t: T) -> Self
     where
         T: RealField,
diff --git a/src/geometry/orthographic.rs b/src/geometry/orthographic.rs
index d6136cec..6cc4dacc 100644
--- a/src/geometry/orthographic.rs
+++ b/src/geometry/orthographic.rs
@@ -188,6 +188,7 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_relative_eq!(proj.as_matrix() * inv, Matrix4::identity());
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse(&self) -> Matrix4<T> {
         let mut res = self.to_homogeneous();
 
@@ -221,6 +222,7 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_eq!(proj.to_homogeneous(), expected);
     /// ```
     #[inline]
+    #[must_use]
     pub fn to_homogeneous(&self) -> Matrix4<T> {
         self.matrix
     }
@@ -240,6 +242,7 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_eq!(*proj.as_matrix(), expected);
     /// ```
     #[inline]
+    #[must_use]
     pub fn as_matrix(&self) -> &Matrix4<T> {
         &self.matrix
     }
@@ -253,6 +256,7 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_eq!(proj.as_projective().to_homogeneous(), proj.to_homogeneous());
     /// ```
     #[inline]
+    #[must_use]
     pub fn as_projective(&self) -> &Projective3<T> {
         unsafe { mem::transmute(self) }
     }
@@ -266,6 +270,7 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_eq!(proj.to_projective().to_homogeneous(), proj.to_homogeneous());
     /// ```
     #[inline]
+    #[must_use]
     pub fn to_projective(&self) -> Projective3<T> {
         Projective3::from_matrix_unchecked(self.matrix)
     }
@@ -310,6 +315,7 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_relative_eq!(proj.left(), 10.0, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn left(&self) -> T {
         (-T::one() - self.matrix[(0, 3)]) / self.matrix[(0, 0)]
     }
@@ -326,6 +332,7 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_relative_eq!(proj.right(), 1.0, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn right(&self) -> T {
         (T::one() - self.matrix[(0, 3)]) / self.matrix[(0, 0)]
     }
@@ -342,6 +349,7 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_relative_eq!(proj.bottom(), 20.0, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn bottom(&self) -> T {
         (-T::one() - self.matrix[(1, 3)]) / self.matrix[(1, 1)]
     }
@@ -358,6 +366,7 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_relative_eq!(proj.top(), 2.0, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn top(&self) -> T {
         (T::one() - self.matrix[(1, 3)]) / self.matrix[(1, 1)]
     }
@@ -374,6 +383,7 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_relative_eq!(proj.znear(), 1000.0, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn znear(&self) -> T {
         (T::one() + self.matrix[(2, 3)]) / self.matrix[(2, 2)]
     }
@@ -390,6 +400,7 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_relative_eq!(proj.zfar(), 0.1, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn zfar(&self) -> T {
         (-T::one() + self.matrix[(2, 3)]) / self.matrix[(2, 2)]
     }
@@ -422,6 +433,7 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_relative_eq!(proj.project_point(&p8), Point3::new( 1.0,  1.0,  1.0));
     /// ```
     #[inline]
+    #[must_use]
     pub fn project_point(&self, p: &Point3<T>) -> Point3<T> {
         Point3::new(
             self.matrix[(0, 0)] * p[0] + self.matrix[(0, 3)],
@@ -457,6 +469,7 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_relative_eq!(proj.unproject_point(&p8), Point3::new(10.0, 20.0, -1000.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn unproject_point(&self, p: &Point3<T>) -> Point3<T> {
         Point3::new(
             (p[0] - self.matrix[(0, 3)]) / self.matrix[(0, 0)],
@@ -485,6 +498,7 @@ impl<T: RealField> Orthographic3<T> {
     /// assert_relative_eq!(proj.project_vector(&v3), Vector3::z() * -2.0 / 999.9);
     /// ```
     #[inline]
+    #[must_use]
     pub fn project_vector<SB>(&self, p: &Vector<T, U3, SB>) -> Vector3<T>
     where
         SB: Storage<T, U3>,
diff --git a/src/geometry/perspective.rs b/src/geometry/perspective.rs
index fac8fa99..65aee5fd 100644
--- a/src/geometry/perspective.rs
+++ b/src/geometry/perspective.rs
@@ -104,6 +104,7 @@ impl<T: RealField> Perspective3<T> {
 
     /// Retrieves the inverse of the underlying homogeneous matrix.
     #[inline]
+    #[must_use]
     pub fn inverse(&self) -> Matrix4<T> {
         let mut res = self.to_homogeneous();
 
@@ -123,24 +124,28 @@ impl<T: RealField> Perspective3<T> {
 
     /// Computes the corresponding homogeneous matrix.
     #[inline]
+    #[must_use]
     pub fn to_homogeneous(&self) -> Matrix4<T> {
         self.matrix.clone_owned()
     }
 
     /// A reference to the underlying homogeneous transformation matrix.
     #[inline]
+    #[must_use]
     pub fn as_matrix(&self) -> &Matrix4<T> {
         &self.matrix
     }
 
     /// A reference to this transformation seen as a `Projective3`.
     #[inline]
+    #[must_use]
     pub fn as_projective(&self) -> &Projective3<T> {
         unsafe { mem::transmute(self) }
     }
 
     /// This transformation seen as a `Projective3`.
     #[inline]
+    #[must_use]
     pub fn to_projective(&self) -> Projective3<T> {
         Projective3::from_matrix_unchecked(self.matrix)
     }
@@ -161,18 +166,21 @@ impl<T: RealField> Perspective3<T> {
 
     /// Gets the `width / height` aspect ratio of the view frustum.
     #[inline]
+    #[must_use]
     pub fn aspect(&self) -> T {
         self.matrix[(1, 1)] / self.matrix[(0, 0)]
     }
 
     /// Gets the y field of view of the view frustum.
     #[inline]
+    #[must_use]
     pub fn fovy(&self) -> T {
         (T::one() / self.matrix[(1, 1)]).atan() * crate::convert(2.0)
     }
 
     /// Gets the near plane offset of the view frustum.
     #[inline]
+    #[must_use]
     pub fn znear(&self) -> T {
         let ratio = (-self.matrix[(2, 2)] + T::one()) / (-self.matrix[(2, 2)] - T::one());
 
@@ -182,6 +190,7 @@ impl<T: RealField> Perspective3<T> {
 
     /// Gets the far plane offset of the view frustum.
     #[inline]
+    #[must_use]
     pub fn zfar(&self) -> T {
         let ratio = (-self.matrix[(2, 2)] + T::one()) / (-self.matrix[(2, 2)] - T::one());
 
@@ -193,6 +202,7 @@ impl<T: RealField> Perspective3<T> {
     // TODO: when we get specialization, specialize the Mul impl instead.
     /// Projects a point. Faster than matrix multiplication.
     #[inline]
+    #[must_use]
     pub fn project_point(&self, p: &Point3<T>) -> Point3<T> {
         let inverse_denom = -T::one() / p[2];
         Point3::new(
@@ -204,6 +214,7 @@ impl<T: RealField> Perspective3<T> {
 
     /// Un-projects a point. Faster than multiplication by the matrix inverse.
     #[inline]
+    #[must_use]
     pub fn unproject_point(&self, p: &Point3<T>) -> Point3<T> {
         let inverse_denom = self.matrix[(2, 3)] / (p[2] + self.matrix[(2, 2)]);
 
@@ -217,6 +228,7 @@ impl<T: RealField> Perspective3<T> {
     // TODO: when we get specialization, specialize the Mul impl instead.
     /// Projects a vector. Faster than matrix multiplication.
     #[inline]
+    #[must_use]
     pub fn project_vector<SB>(&self, p: &Vector<T, U3, SB>) -> Vector3<T>
     where
         SB: Storage<T, U3>,
diff --git a/src/geometry/point.rs b/src/geometry/point.rs
index 60fa9156..d25b4d7c 100644
--- a/src/geometry/point.rs
+++ b/src/geometry/point.rs
@@ -122,6 +122,7 @@ impl<T: Scalar, const D: usize> Point<T, D> {
     /// assert_eq!(p.map(|e| e as u32), Point3::new(1, 2, 3));
     /// ```
     #[inline]
+    #[must_use]
     pub fn map<T2: Scalar, F: FnMut(T) -> T2>(&self, f: F) -> Point<T2, D> {
         self.coords.map(f).into()
     }
@@ -161,6 +162,7 @@ impl<T: Scalar, const D: usize> Point<T, D> {
     /// assert_eq!(p.to_homogeneous(), Vector4::new(10.0, 20.0, 30.0, 1.0));
     /// ```
     #[inline]
+    #[must_use]
     pub fn to_homogeneous(&self) -> OVector<T, DimNameSum<Const<D>, U1>>
     where
         T: One,
@@ -199,6 +201,7 @@ impl<T: Scalar, const D: usize> Point<T, D> {
     /// assert_eq!(p.len(), 3);
     /// ```
     #[inline]
+    #[must_use]
     pub fn len(&self) -> usize {
         self.coords.len()
     }
@@ -212,6 +215,7 @@ impl<T: Scalar, const D: usize> Point<T, D> {
     /// assert!(!p.is_empty());
     /// ```
     #[inline]
+    #[must_use]
     pub fn is_empty(&self) -> bool {
         self.len() == 0
     }
@@ -246,6 +250,7 @@ impl<T: Scalar, const D: usize> Point<T, D> {
 
     /// Gets a reference to i-th element of this point without bound-checking.
     #[inline]
+    #[must_use]
     pub unsafe fn get_unchecked(&self, i: usize) -> &T {
         self.coords.vget_unchecked(i)
     }
@@ -272,6 +277,7 @@ impl<T: Scalar, const D: usize> Point<T, D> {
 
     /// Gets a mutable reference to i-th element of this point without bound-checking.
     #[inline]
+    #[must_use]
     pub unsafe fn get_unchecked_mut(&mut self, i: usize) -> &mut T {
         self.coords.vget_unchecked_mut(i)
     }
@@ -378,18 +384,21 @@ impl<T: Scalar + PartialOrd, const D: usize> PartialOrd for Point<T, D> {
 impl<T: Scalar + SimdPartialOrd, const D: usize> Point<T, D> {
     /// Computes the infimum (aka. componentwise min) of two points.
     #[inline]
+    #[must_use]
     pub fn inf(&self, other: &Self) -> Point<T, D> {
         self.coords.inf(&other.coords).into()
     }
 
     /// Computes the supremum (aka. componentwise max) of two points.
     #[inline]
+    #[must_use]
     pub fn sup(&self, other: &Self) -> Point<T, D> {
         self.coords.sup(&other.coords).into()
     }
 
     /// Computes the (infimum, supremum) of two points.
     #[inline]
+    #[must_use]
     pub fn inf_sup(&self, other: &Self) -> (Point<T, D>, Point<T, D>) {
         let (inf, sup) = self.coords.inf_sup(&other.coords);
         (inf.into(), sup.into())
diff --git a/src/geometry/quaternion.rs b/src/geometry/quaternion.rs
index 6852335b..79fda90a 100755
--- a/src/geometry/quaternion.rs
+++ b/src/geometry/quaternion.rs
@@ -191,6 +191,7 @@ where
 
     /// The imaginary part of this quaternion.
     #[inline]
+    #[must_use]
     pub fn imag(&self) -> Vector3<T> {
         self.coords.xyz()
     }
@@ -223,6 +224,7 @@ where
     /// assert_eq!(q1.lerp(&q2, 0.1), Quaternion::new(1.9, 3.8, 5.7, 7.6));
     /// ```
     #[inline]
+    #[must_use]
     pub fn lerp(&self, other: &Self, t: T) -> Self {
         self * (T::one() - t) + other * t
     }
@@ -238,6 +240,7 @@ where
     /// assert_eq!(q.vector()[2], 4.0);
     /// ```
     #[inline]
+    #[must_use]
     pub fn vector(&self) -> MatrixSlice<T, U3, U1, RStride<T, U4, U1>, CStride<T, U4, U1>> {
         self.coords.fixed_rows::<3>(0)
     }
@@ -251,6 +254,7 @@ where
     /// assert_eq!(q.scalar(), 1.0);
     /// ```
     #[inline]
+    #[must_use]
     pub fn scalar(&self) -> T {
         self.coords[3]
     }
@@ -266,6 +270,7 @@ where
     /// assert_eq!(*q.as_vector(), Vector4::new(2.0, 3.0, 4.0, 1.0));
     /// ```
     #[inline]
+    #[must_use]
     pub fn as_vector(&self) -> &Vector4<T> {
         &self.coords
     }
@@ -280,6 +285,7 @@ where
     /// assert_relative_eq!(q.norm(), 5.47722557, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn norm(&self) -> T {
         self.coords.norm()
     }
@@ -297,6 +303,7 @@ where
     /// assert_relative_eq!(q.magnitude(), 5.47722557, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn magnitude(&self) -> T {
         self.norm()
     }
@@ -310,6 +317,7 @@ where
     /// assert_eq!(q.magnitude_squared(), 30.0);
     /// ```
     #[inline]
+    #[must_use]
     pub fn norm_squared(&self) -> T {
         self.coords.norm_squared()
     }
@@ -326,6 +334,7 @@ where
     /// assert_eq!(q.magnitude_squared(), 30.0);
     /// ```
     #[inline]
+    #[must_use]
     pub fn magnitude_squared(&self) -> T {
         self.norm_squared()
     }
@@ -340,6 +349,7 @@ where
     /// assert_eq!(q1.dot(&q2), 70.0);
     /// ```
     #[inline]
+    #[must_use]
     pub fn dot(&self, rhs: &Self) -> T {
         self.coords.dot(&rhs.coords)
     }
@@ -409,6 +419,7 @@ where
     /// let result = a.inner(&b);
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-5);
     #[inline]
+    #[must_use]
     pub fn inner(&self, other: &Self) -> Self {
         (self * other + other * self).half()
     }
@@ -428,6 +439,7 @@ where
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-5);
     /// ```
     #[inline]
+    #[must_use]
     pub fn outer(&self, other: &Self) -> Self {
         #[allow(clippy::eq_op)]
         (self * other - other * self).half()
@@ -448,6 +460,7 @@ where
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-5);
     /// ```
     #[inline]
+    #[must_use]
     pub fn project(&self, other: &Self) -> Option<Self>
     where
         T: RealField,
@@ -470,6 +483,7 @@ where
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-5);
     /// ```
     #[inline]
+    #[must_use]
     pub fn reject(&self, other: &Self) -> Option<Self>
     where
         T: RealField,
@@ -492,6 +506,7 @@ where
     /// assert_eq!(half_ang, f32::consts::FRAC_PI_2);
     /// assert_eq!(axis, Some(Vector3::x_axis()));
     /// ```
+    #[must_use]
     pub fn polar_decomposition(&self) -> (T, T, Option<Unit<Vector3<T>>>)
     where
         T: RealField,
@@ -519,6 +534,7 @@ where
     /// assert_relative_eq!(q.ln(), Quaternion::new(1.683647, 1.190289, 0.0, 0.0), epsilon = 1.0e-6)
     /// ```
     #[inline]
+    #[must_use]
     pub fn ln(&self) -> Self {
         let n = self.norm();
         let v = self.vector();
@@ -537,6 +553,7 @@ where
     /// assert_relative_eq!(q.exp(), Quaternion::new(2.0, 5.0, 0.0, 0.0), epsilon = 1.0e-5)
     /// ```
     #[inline]
+    #[must_use]
     pub fn exp(&self) -> Self {
         self.exp_eps(T::simd_default_epsilon())
     }
@@ -556,6 +573,7 @@ where
     /// assert_eq!(q.exp_eps(1.0e-6), Quaternion::identity());
     /// ```
     #[inline]
+    #[must_use]
     pub fn exp_eps(&self, eps: T) -> Self {
         let v = self.vector();
         let nn = v.norm_squared();
@@ -579,6 +597,7 @@ where
     /// assert_relative_eq!(q.powf(1.5), Quaternion::new( -6.2576659, 4.1549037, 6.2323556, 8.3098075), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn powf(&self, n: T) -> Self {
         (self.ln() * n).exp()
     }
@@ -674,18 +693,21 @@ where
 
     /// Calculates square of a quaternion.
     #[inline]
+    #[must_use]
     pub fn squared(&self) -> Self {
         self * self
     }
 
     /// Divides quaternion into two.
     #[inline]
+    #[must_use]
     pub fn half(&self) -> Self {
         self / crate::convert(2.0f64)
     }
 
     /// Calculates square root.
     #[inline]
+    #[must_use]
     pub fn sqrt(&self) -> Self {
         self.powf(crate::convert(0.5))
     }
@@ -694,12 +716,14 @@ where
     ///
     /// A quaternion is pure if it has no real part (`self.w == 0.0`).
     #[inline]
+    #[must_use]
     pub fn is_pure(&self) -> bool {
         self.w.is_zero()
     }
 
     /// Convert quaternion to pure quaternion.
     #[inline]
+    #[must_use]
     pub fn pure(&self) -> Self {
         Self::from_imag(self.imag())
     }
@@ -708,6 +732,7 @@ where
     ///
     /// Calculates B<sup>-1</sup> * A where A = self, B = other.
     #[inline]
+    #[must_use]
     pub fn left_div(&self, other: &Self) -> Option<Self>
     where
         T: RealField,
@@ -730,6 +755,7 @@ where
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-7);
     /// ```
     #[inline]
+    #[must_use]
     pub fn right_div(&self, other: &Self) -> Option<Self>
     where
         T: RealField,
@@ -749,6 +775,7 @@ where
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-7);
     /// ```
     #[inline]
+    #[must_use]
     pub fn cos(&self) -> Self {
         let z = self.imag().magnitude();
         let w = -self.w.simd_sin() * z.simd_sinhc();
@@ -766,6 +793,7 @@ where
     /// assert_relative_eq!(input, result, epsilon = 1.0e-7);
     /// ```
     #[inline]
+    #[must_use]
     pub fn acos(&self) -> Self {
         let u = Self::from_imag(self.imag().normalize());
         let identity = Self::identity();
@@ -787,6 +815,7 @@ where
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-7);
     /// ```
     #[inline]
+    #[must_use]
     pub fn sin(&self) -> Self {
         let z = self.imag().magnitude();
         let w = self.w.simd_cos() * z.simd_sinhc();
@@ -804,6 +833,7 @@ where
     /// assert_relative_eq!(input, result, epsilon = 1.0e-7);
     /// ```
     #[inline]
+    #[must_use]
     pub fn asin(&self) -> Self {
         let u = Self::from_imag(self.imag().normalize());
         let identity = Self::identity();
@@ -825,6 +855,7 @@ where
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-7);
     /// ```
     #[inline]
+    #[must_use]
     pub fn tan(&self) -> Self
     where
         T: RealField,
@@ -843,6 +874,7 @@ where
     /// assert_relative_eq!(input, result, epsilon = 1.0e-7);
     /// ```
     #[inline]
+    #[must_use]
     pub fn atan(&self) -> Self
     where
         T: RealField,
@@ -867,6 +899,7 @@ where
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-7);
     /// ```
     #[inline]
+    #[must_use]
     pub fn sinh(&self) -> Self {
         (self.exp() - (-self).exp()).half()
     }
@@ -883,6 +916,7 @@ where
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-7);
     /// ```
     #[inline]
+    #[must_use]
     pub fn asinh(&self) -> Self {
         let identity = Self::identity();
         (self + (identity + self.squared()).sqrt()).ln()
@@ -900,6 +934,7 @@ where
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-7);
     /// ```
     #[inline]
+    #[must_use]
     pub fn cosh(&self) -> Self {
         (self.exp() + (-self).exp()).half()
     }
@@ -916,6 +951,7 @@ where
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-7);
     /// ```
     #[inline]
+    #[must_use]
     pub fn acosh(&self) -> Self {
         let identity = Self::identity();
         (self + (self + identity).sqrt() * (self - identity).sqrt()).ln()
@@ -933,6 +969,7 @@ where
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-7);
     /// ```
     #[inline]
+    #[must_use]
     pub fn tanh(&self) -> Self
     where
         T: RealField,
@@ -952,6 +989,7 @@ where
     /// assert_relative_eq!(expected, result, epsilon = 1.0e-7);
     /// ```
     #[inline]
+    #[must_use]
     pub fn atanh(&self) -> Self {
         let identity = Self::identity();
         ((identity + self).ln() - (identity - self).ln()).half()
@@ -1069,6 +1107,7 @@ where
     /// assert_eq!(rot.angle(), 1.78);
     /// ```
     #[inline]
+    #[must_use]
     pub fn angle(&self) -> T {
         let w = self.quaternion().scalar().simd_abs();
         self.quaternion().imag().norm().simd_atan2(w) * crate::convert(2.0f64)
@@ -1085,6 +1124,7 @@ where
     /// assert_eq!(*axis.quaternion(), Quaternion::new(1.0, 0.0, 0.0, 0.0));
     /// ```
     #[inline]
+    #[must_use]
     pub fn quaternion(&self) -> &Quaternion<T> {
         self.as_ref()
     }
@@ -1133,6 +1173,7 @@ where
     /// assert_relative_eq!(rot1.angle_to(&rot2), 1.0045657, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn angle_to(&self, other: &Self) -> T {
         let delta = self.rotation_to(other);
         delta.angle()
@@ -1152,6 +1193,7 @@ where
     /// assert_relative_eq!(rot_to * rot1, rot2, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn rotation_to(&self, other: &Self) -> Self {
         other / self
     }
@@ -1168,6 +1210,7 @@ where
     /// assert_eq!(q1.lerp(&q2, 0.1), Quaternion::new(0.9, 0.1, 0.0, 0.0));
     /// ```
     #[inline]
+    #[must_use]
     pub fn lerp(&self, other: &Self, t: T) -> Quaternion<T> {
         self.as_ref().lerp(other.as_ref(), t)
     }
@@ -1184,6 +1227,7 @@ where
     /// assert_eq!(q1.nlerp(&q2, 0.1), UnitQuaternion::new_normalize(Quaternion::new(0.9, 0.1, 0.0, 0.0)));
     /// ```
     #[inline]
+    #[must_use]
     pub fn nlerp(&self, other: &Self, t: T) -> Self {
         let mut res = self.lerp(other, t);
         let _ = res.normalize_mut();
@@ -1209,6 +1253,7 @@ where
     /// assert_eq!(q.euler_angles(), (std::f32::consts::FRAC_PI_2, 0.0, 0.0));
     /// ```
     #[inline]
+    #[must_use]
     pub fn slerp(&self, other: &Self, t: T) -> Self
     where
         T: RealField,
@@ -1228,6 +1273,7 @@ where
     /// * `epsilon`: the value below which the sinus of the angle separating both quaternion
     /// must be to return `None`.
     #[inline]
+    #[must_use]
     pub fn try_slerp(&self, other: &Self, t: T, epsilon: T) -> Option<Self>
     where
         T: RealField,
@@ -1287,6 +1333,7 @@ where
     /// assert!(rot.axis().is_none());
     /// ```
     #[inline]
+    #[must_use]
     pub fn axis(&self) -> Option<Unit<Vector3<T>>>
     where
         T: RealField,
@@ -1311,6 +1358,7 @@ where
     /// assert_relative_eq!(rot.scaled_axis(), axisangle, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn scaled_axis(&self) -> Vector3<T>
     where
         T: RealField,
@@ -1339,6 +1387,7 @@ where
     /// assert!(rot.axis_angle().is_none());
     /// ```
     #[inline]
+    #[must_use]
     pub fn axis_angle(&self) -> Option<(Unit<Vector3<T>>, T)>
     where
         T: RealField,
@@ -1350,6 +1399,7 @@ where
     ///
     /// Note that this function yields a `Quaternion<T>` because it loses the unit property.
     #[inline]
+    #[must_use]
     pub fn exp(&self) -> Quaternion<T> {
         self.as_ref().exp()
     }
@@ -1369,6 +1419,7 @@ where
     /// assert_relative_eq!(q.ln().vector().into_owned(), axisangle, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn ln(&self) -> Quaternion<T>
     where
         T: RealField,
@@ -1397,6 +1448,7 @@ where
     /// assert_eq!(pow.angle(), 2.4);
     /// ```
     #[inline]
+    #[must_use]
     pub fn powf(&self, n: T) -> Self
     where
         T: RealField,
@@ -1425,6 +1477,7 @@ where
     /// assert_relative_eq!(*rot.matrix(), expected, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn to_rotation_matrix(&self) -> Rotation<T, 3> {
         let i = self.as_ref()[0];
         let j = self.as_ref()[1];
@@ -1482,6 +1535,7 @@ where
     /// assert_relative_eq!(euler.2, 0.3, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn euler_angles(&self) -> (T, T, T)
     where
         T: RealField,
@@ -1506,6 +1560,7 @@ where
     /// assert_relative_eq!(rot.to_homogeneous(), expected, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn to_homogeneous(&self) -> Matrix4<T> {
         self.to_rotation_matrix().to_homogeneous()
     }
@@ -1526,6 +1581,7 @@ where
     /// assert_relative_eq!(transformed_point, Point3::new(3.0, 2.0, -1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn transform_point(&self, pt: &Point3<T>) -> Point3<T> {
         self * pt
     }
@@ -1546,6 +1602,7 @@ where
     /// assert_relative_eq!(transformed_vector, Vector3::new(3.0, 2.0, -1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn transform_vector(&self, v: &Vector3<T>) -> Vector3<T> {
         self * v
     }
@@ -1566,6 +1623,7 @@ where
     /// assert_relative_eq!(transformed_point, Point3::new(-3.0, 2.0, 1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_point(&self, pt: &Point3<T>) -> Point3<T> {
         // TODO: would it be useful performancewise not to call inverse explicitly (i-e. implement
         // the inverse transformation explicitly here) ?
@@ -1588,6 +1646,7 @@ where
     /// assert_relative_eq!(transformed_vector, Vector3::new(-3.0, 2.0, 1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_vector(&self, v: &Vector3<T>) -> Vector3<T> {
         self.inverse() * v
     }
@@ -1608,6 +1667,7 @@ where
     /// assert_relative_eq!(transformed_vector, -Vector3::y_axis(), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_unit_vector(&self, v: &Unit<Vector3<T>>) -> Unit<Vector3<T>> {
         self.inverse() * v
     }
@@ -1616,6 +1676,7 @@ where
     ///
     /// This is faster, but approximate, way to compute `UnitQuaternion::new(axisangle) * self`.
     #[inline]
+    #[must_use]
     pub fn append_axisangle_linearized(&self, axisangle: &Vector3<T>) -> Self {
         let half: T = crate::convert(0.5);
         let q1 = self.into_inner();
diff --git a/src/geometry/reflection.rs b/src/geometry/reflection.rs
index e48c700a..dcc754c2 100644
--- a/src/geometry/reflection.rs
+++ b/src/geometry/reflection.rs
@@ -34,6 +34,7 @@ impl<T: ComplexField, D: Dim, S: Storage<T, D>> Reflection<T, D, S> {
     }
 
     /// The reflexion axis.
+    #[must_use]
     pub fn axis(&self) -> &Vector<T, D, S> {
         &self.axis
     }
diff --git a/src/geometry/rotation.rs b/src/geometry/rotation.rs
index 7b4d609a..f3127fb9 100755
--- a/src/geometry/rotation.rs
+++ b/src/geometry/rotation.rs
@@ -185,6 +185,7 @@ impl<T: Scalar, const D: usize> Rotation<T, D> {
     /// assert_eq!(*rot.matrix(), expected);
     /// ```
     #[inline]
+    #[must_use]
     pub fn matrix(&self) -> &SMatrix<T, D, D> {
         &self.matrix
     }
@@ -262,6 +263,7 @@ impl<T: Scalar, const D: usize> Rotation<T, D> {
     /// assert_eq!(rot.to_homogeneous(), expected);
     /// ```
     #[inline]
+    #[must_use]
     pub fn to_homogeneous(&self) -> OMatrix<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>
     where
         T: Zero + One,
@@ -403,6 +405,7 @@ where
     /// assert_relative_eq!(transformed_point, Point3::new(3.0, 2.0, -1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn transform_point(&self, pt: &Point<T, D>) -> Point<T, D> {
         self * pt
     }
@@ -422,6 +425,7 @@ where
     /// assert_relative_eq!(transformed_vector, Vector3::new(3.0, 2.0, -1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn transform_vector(&self, v: &SVector<T, D>) -> SVector<T, D> {
         self * v
     }
@@ -441,6 +445,7 @@ where
     /// assert_relative_eq!(transformed_point, Point3::new(-3.0, 2.0, 1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_point(&self, pt: &Point<T, D>) -> Point<T, D> {
         Point::from(self.inverse_transform_vector(&pt.coords))
     }
@@ -460,6 +465,7 @@ where
     /// assert_relative_eq!(transformed_vector, Vector3::new(-3.0, 2.0, 1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_vector(&self, v: &SVector<T, D>) -> SVector<T, D> {
         self.matrix().tr_mul(v)
     }
@@ -479,6 +485,7 @@ where
     /// assert_relative_eq!(transformed_vector, -Vector3::y_axis(), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_unit_vector(&self, v: &Unit<SVector<T, D>>) -> Unit<SVector<T, D>> {
         Unit::new_unchecked(self.inverse_transform_vector(&**v))
     }
diff --git a/src/geometry/rotation_interpolation.rs b/src/geometry/rotation_interpolation.rs
index bf1c9094..6b5a5c11 100644
--- a/src/geometry/rotation_interpolation.rs
+++ b/src/geometry/rotation_interpolation.rs
@@ -18,6 +18,7 @@ impl<T: SimdRealField> Rotation2<T> {
     /// assert_relative_eq!(rot.angle(), std::f32::consts::FRAC_PI_2);
     /// ```
     #[inline]
+    #[must_use]
     pub fn slerp(&self, other: &Self, t: T) -> Self
     where
         T::Element: SimdRealField,
@@ -47,6 +48,7 @@ impl<T: SimdRealField> Rotation3<T> {
     /// assert_eq!(q.euler_angles(), (std::f32::consts::FRAC_PI_2, 0.0, 0.0));
     /// ```
     #[inline]
+    #[must_use]
     pub fn slerp(&self, other: &Self, t: T) -> Self
     where
         T: RealField,
@@ -67,6 +69,7 @@ impl<T: SimdRealField> Rotation3<T> {
     /// * `epsilon`: the value below which the sinus of the angle separating both rotations
     /// must be to return `None`.
     #[inline]
+    #[must_use]
     pub fn try_slerp(&self, other: &Self, t: T, epsilon: T) -> Option<Self>
     where
         T: RealField,
diff --git a/src/geometry/rotation_specialization.rs b/src/geometry/rotation_specialization.rs
index bb71f4f8..a89f5845 100644
--- a/src/geometry/rotation_specialization.rs
+++ b/src/geometry/rotation_specialization.rs
@@ -186,6 +186,7 @@ impl<T: SimdRealField> Rotation2<T> {
     /// assert_relative_eq!(rot_to.inverse() * rot2, rot1);
     /// ```
     #[inline]
+    #[must_use]
     pub fn rotation_to(&self, other: &Self) -> Self {
         other * self.inverse()
     }
@@ -215,6 +216,7 @@ impl<T: SimdRealField> Rotation2<T> {
     /// assert_relative_eq!(pow.angle(), 2.0 * 0.78);
     /// ```
     #[inline]
+    #[must_use]
     pub fn powf(&self, n: T) -> Self {
         Self::new(self.angle() * n)
     }
@@ -232,6 +234,7 @@ impl<T: SimdRealField> Rotation2<T> {
     /// assert_relative_eq!(rot.angle(), 1.78);
     /// ```
     #[inline]
+    #[must_use]
     pub fn angle(&self) -> T {
         self.matrix()[(1, 0)].simd_atan2(self.matrix()[(0, 0)])
     }
@@ -247,6 +250,7 @@ impl<T: SimdRealField> Rotation2<T> {
     /// assert_relative_eq!(rot1.angle_to(&rot2), 1.6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn angle_to(&self, other: &Self) -> T {
         self.rotation_to(other).angle()
     }
@@ -256,6 +260,7 @@ impl<T: SimdRealField> Rotation2<T> {
     /// This is generally used in the context of generic programming. Using
     /// the `.angle()` method instead is more common.
     #[inline]
+    #[must_use]
     pub fn scaled_axis(&self) -> SVector<T, 1> {
         Vector1::new(self.angle())
     }
@@ -640,6 +645,7 @@ where
     /// assert_relative_eq!(rot_to * rot1, rot2, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn rotation_to(&self, other: &Self) -> Self {
         other * self.inverse()
     }
@@ -659,6 +665,7 @@ where
     /// assert_eq!(pow.angle(), 2.4);
     /// ```
     #[inline]
+    #[must_use]
     pub fn powf(&self, n: T) -> Self
     where
         T: RealField,
@@ -765,6 +772,7 @@ impl<T: SimdRealField> Rotation3<T> {
     /// assert_relative_eq!(rot.angle(), 1.78);
     /// ```
     #[inline]
+    #[must_use]
     pub fn angle(&self) -> T {
         ((self.matrix()[(0, 0)] + self.matrix()[(1, 1)] + self.matrix()[(2, 2)] - T::one())
             / crate::convert(2.0))
@@ -787,6 +795,7 @@ impl<T: SimdRealField> Rotation3<T> {
     /// assert!(rot.axis().is_none());
     /// ```
     #[inline]
+    #[must_use]
     pub fn axis(&self) -> Option<Unit<Vector3<T>>>
     where
         T: RealField,
@@ -811,6 +820,7 @@ impl<T: SimdRealField> Rotation3<T> {
     /// assert_relative_eq!(rot.scaled_axis(), axisangle, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn scaled_axis(&self) -> Vector3<T>
     where
         T: RealField,
@@ -842,6 +852,7 @@ impl<T: SimdRealField> Rotation3<T> {
     /// assert!(rot.axis_angle().is_none());
     /// ```
     #[inline]
+    #[must_use]
     pub fn axis_angle(&self) -> Option<(Unit<Vector3<T>>, T)>
     where
         T: RealField,
@@ -864,6 +875,7 @@ impl<T: SimdRealField> Rotation3<T> {
     /// assert_relative_eq!(rot1.angle_to(&rot2), 1.0045657, epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn angle_to(&self, other: &Self) -> T
     where
         T::Element: SimdRealField,
@@ -896,6 +908,7 @@ impl<T: SimdRealField> Rotation3<T> {
     /// assert_relative_eq!(euler.1, 0.2, epsilon = 1.0e-6);
     /// assert_relative_eq!(euler.2, 0.3, epsilon = 1.0e-6);
     /// ```
+    #[must_use]
     pub fn euler_angles(&self) -> (T, T, T)
     where
         T: RealField,
diff --git a/src/geometry/similarity.rs b/src/geometry/similarity.rs
index b6e2c379..516b95c5 100755
--- a/src/geometry/similarity.rs
+++ b/src/geometry/similarity.rs
@@ -122,6 +122,7 @@ where
 impl<T: Scalar, R, const D: usize> Similarity<T, R, D> {
     /// The scaling factor of this similarity transformation.
     #[inline]
+    #[must_use]
     pub fn scaling(&self) -> T {
         self.scaling.inlined_clone()
     }
@@ -248,6 +249,7 @@ where
     /// assert_relative_eq!(transformed_point, Point3::new(19.0, 17.0, -9.0), epsilon = 1.0e-5);
     /// ```
     #[inline]
+    #[must_use]
     pub fn transform_point(&self, pt: &Point<T, D>) -> Point<T, D> {
         self * pt
     }
@@ -269,6 +271,7 @@ where
     /// assert_relative_eq!(transformed_vector, Vector3::new(18.0, 15.0, -12.0), epsilon = 1.0e-5);
     /// ```
     #[inline]
+    #[must_use]
     pub fn transform_vector(&self, v: &SVector<T, D>) -> SVector<T, D> {
         self * v
     }
@@ -289,6 +292,7 @@ where
     /// assert_relative_eq!(transformed_point, Point3::new(-1.5, 1.5, 1.5), epsilon = 1.0e-5);
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_point(&self, pt: &Point<T, D>) -> Point<T, D> {
         self.isometry.inverse_transform_point(pt) / self.scaling()
     }
@@ -309,6 +313,7 @@ where
     /// assert_relative_eq!(transformed_vector, Vector3::new(-3.0, 2.5, 2.0), epsilon = 1.0e-5);
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_vector(&self, v: &SVector<T, D>) -> SVector<T, D> {
         self.isometry.inverse_transform_vector(v) / self.scaling()
     }
@@ -321,6 +326,7 @@ where
 impl<T: SimdRealField, R, const D: usize> Similarity<T, R, D> {
     /// Converts this similarity into its equivalent homogeneous transformation matrix.
     #[inline]
+    #[must_use]
     pub fn to_homogeneous(&self) -> OMatrix<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>
     where
         Const<D>: DimNameAdd<U1>,
diff --git a/src/geometry/swizzle.rs b/src/geometry/swizzle.rs
index 6e4c6df8..0ad51f00 100644
--- a/src/geometry/swizzle.rs
+++ b/src/geometry/swizzle.rs
@@ -8,6 +8,7 @@ macro_rules! impl_swizzle {
             $(
                 /// Builds a new point from components of `self`.
                 #[inline]
+                #[must_use]
                 pub fn $name(&self) -> $Result<T>
                  where <Const<D> as ToTypenum>::Typenum: Cmp<typenum::$BaseDim, Output=Greater> {
                     $Result::new($(self[$i].inlined_clone()),*)
diff --git a/src/geometry/transform.rs b/src/geometry/transform.rs
index c7fce902..8c294b15 100755
--- a/src/geometry/transform.rs
+++ b/src/geometry/transform.rs
@@ -301,6 +301,7 @@ where
     /// assert_eq!(*t.matrix(), m);
     /// ```
     #[inline]
+    #[must_use]
     pub fn matrix(&self) -> &OMatrix<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>> {
         &self.matrix
     }
@@ -367,6 +368,7 @@ where
     /// assert_eq!(t.into_inner(), m);
     /// ```
     #[inline]
+    #[must_use]
     pub fn to_homogeneous(&self) -> OMatrix<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>> {
         self.matrix().clone_owned()
     }
@@ -498,6 +500,7 @@ where
     ///
     /// This is the same as the multiplication `self * pt`.
     #[inline]
+    #[must_use]
     pub fn transform_point(&self, pt: &Point<T, D>) -> Point<T, D> {
         self * pt
     }
@@ -507,6 +510,7 @@ where
     ///
     /// This is the same as the multiplication `self * v`.
     #[inline]
+    #[must_use]
     pub fn transform_vector(&self, v: &SVector<T, D>) -> SVector<T, D> {
         self * v
     }
@@ -524,6 +528,7 @@ where
     /// This may be cheaper than inverting the transformation and transforming
     /// the point.
     #[inline]
+    #[must_use]
     pub fn inverse_transform_point(&self, pt: &Point<T, D>) -> Point<T, D> {
         self.clone().inverse() * pt
     }
@@ -532,6 +537,7 @@ where
     /// This may be cheaper than inverting the transformation and transforming
     /// the vector.
     #[inline]
+    #[must_use]
     pub fn inverse_transform_vector(&self, v: &SVector<T, D>) -> SVector<T, D> {
         self.clone().inverse() * v
     }
diff --git a/src/geometry/translation.rs b/src/geometry/translation.rs
index dcb7b603..1476e75b 100755
--- a/src/geometry/translation.rs
+++ b/src/geometry/translation.rs
@@ -190,6 +190,7 @@ impl<T: Scalar, const D: usize> Translation<T, D> {
     /// assert_eq!(t.to_homogeneous(), expected);
     /// ```
     #[inline]
+    #[must_use]
     pub fn to_homogeneous(&self) -> OMatrix<T, DimNameSum<Const<D>, U1>, DimNameSum<Const<D>, U1>>
     where
         T: Zero + One,
@@ -241,6 +242,7 @@ impl<T: Scalar + ClosedAdd, const D: usize> Translation<T, D> {
     /// let transformed_point = t.transform_point(&Point3::new(4.0, 5.0, 6.0));
     /// assert_eq!(transformed_point, Point3::new(5.0, 7.0, 9.0));
     #[inline]
+    #[must_use]
     pub fn transform_point(&self, pt: &Point<T, D>) -> Point<T, D> {
         pt + &self.vector
     }
@@ -256,6 +258,7 @@ impl<T: Scalar + ClosedSub, const D: usize> Translation<T, D> {
     /// let transformed_point = t.inverse_transform_point(&Point3::new(4.0, 5.0, 6.0));
     /// assert_eq!(transformed_point, Point3::new(3.0, 3.0, 3.0));
     #[inline]
+    #[must_use]
     pub fn inverse_transform_point(&self, pt: &Point<T, D>) -> Point<T, D> {
         pt - &self.vector
     }
diff --git a/src/geometry/unit_complex.rs b/src/geometry/unit_complex.rs
index fe2e87ba..be0c8980 100755
--- a/src/geometry/unit_complex.rs
+++ b/src/geometry/unit_complex.rs
@@ -84,6 +84,7 @@ where
     /// assert_eq!(rot.angle(), 1.78);
     /// ```
     #[inline]
+    #[must_use]
     pub fn angle(&self) -> T {
         self.im.simd_atan2(self.re)
     }
@@ -98,6 +99,7 @@ where
     /// assert_eq!(rot.sin_angle(), angle.sin());
     /// ```
     #[inline]
+    #[must_use]
     pub fn sin_angle(&self) -> T {
         self.im
     }
@@ -112,6 +114,7 @@ where
     /// assert_eq!(rot.cos_angle(),angle.cos());
     /// ```
     #[inline]
+    #[must_use]
     pub fn cos_angle(&self) -> T {
         self.re
     }
@@ -121,6 +124,7 @@ where
     /// This is generally used in the context of generic programming. Using
     /// the `.angle()` method instead is more common.
     #[inline]
+    #[must_use]
     pub fn scaled_axis(&self) -> Vector1<T> {
         Vector1::new(self.angle())
     }
@@ -131,6 +135,7 @@ where
     /// the `.angle()` method instead is more common.
     /// Returns `None` if the angle is zero.
     #[inline]
+    #[must_use]
     pub fn axis_angle(&self) -> Option<(Unit<Vector1<T>>, T)>
     where
         T: RealField,
@@ -157,6 +162,7 @@ where
     /// assert_relative_eq!(rot1.angle_to(&rot2), 1.6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn angle_to(&self, other: &Self) -> T {
         let delta = self.rotation_to(other);
         delta.angle()
@@ -254,6 +260,7 @@ where
     /// assert_eq!(rot.to_rotation_matrix(), expected);
     /// ```
     #[inline]
+    #[must_use]
     pub fn to_rotation_matrix(&self) -> Rotation2<T> {
         let r = self.re;
         let i = self.im;
@@ -274,6 +281,7 @@ where
     /// assert_eq!(rot.to_homogeneous(), expected);
     /// ```
     #[inline]
+    #[must_use]
     pub fn to_homogeneous(&self) -> Matrix3<T> {
         self.to_rotation_matrix().to_homogeneous()
     }
@@ -298,6 +306,7 @@ where
     /// assert_relative_eq!(transformed_point, Point2::new(-2.0, 1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn transform_point(&self, pt: &Point2<T>) -> Point2<T> {
         self * pt
     }
@@ -316,6 +325,7 @@ where
     /// assert_relative_eq!(transformed_vector, Vector2::new(-2.0, 1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn transform_vector(&self, v: &Vector2<T>) -> Vector2<T> {
         self * v
     }
@@ -332,6 +342,7 @@ where
     /// assert_relative_eq!(transformed_point, Point2::new(2.0, -1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_point(&self, pt: &Point2<T>) -> Point2<T> {
         // TODO: would it be useful performancewise not to call inverse explicitly (i-e. implement
         // the inverse transformation explicitly here) ?
@@ -350,6 +361,7 @@ where
     /// assert_relative_eq!(transformed_vector, Vector2::new(2.0, -1.0), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_vector(&self, v: &Vector2<T>) -> Vector2<T> {
         self.inverse() * v
     }
@@ -366,6 +378,7 @@ where
     /// assert_relative_eq!(transformed_vector, -Vector2::y_axis(), epsilon = 1.0e-6);
     /// ```
     #[inline]
+    #[must_use]
     pub fn inverse_transform_unit_vector(&self, v: &Unit<Vector2<T>>) -> Unit<Vector2<T>> {
         self.inverse() * v
     }
@@ -392,6 +405,7 @@ where
     /// assert_relative_eq!(rot.angle(), std::f32::consts::FRAC_PI_2);
     /// ```
     #[inline]
+    #[must_use]
     pub fn slerp(&self, other: &Self, t: T) -> Self {
         Self::new(self.angle() * (T::one() - t) + other.angle() * t)
     }
diff --git a/src/geometry/unit_complex_construction.rs b/src/geometry/unit_complex_construction.rs
index 62a9b416..16e605f2 100644
--- a/src/geometry/unit_complex_construction.rs
+++ b/src/geometry/unit_complex_construction.rs
@@ -148,6 +148,7 @@ where
     /// assert_eq!(*rot.complex(), Complex::new(angle.cos(), angle.sin()));
     /// ```
     #[inline]
+    #[must_use]
     pub fn complex(&self) -> &Complex<T> {
         self.as_ref()
     }
@@ -244,6 +245,7 @@ where
     /// assert_relative_eq!(rot_to.inverse() * rot2, rot1);
     /// ```
     #[inline]
+    #[must_use]
     pub fn rotation_to(&self, other: &Self) -> Self {
         other / self
     }
@@ -262,6 +264,7 @@ where
     /// assert_relative_eq!(pow.angle(), 2.0 * 0.78);
     /// ```
     #[inline]
+    #[must_use]
     pub fn powf(&self, n: T) -> Self {
         Self::from_angle(self.angle() * n)
     }
diff --git a/src/lib.rs b/src/lib.rs
index 8378b272..a978d5f5 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -184,6 +184,7 @@ pub fn zero<T: Zero>() -> T {
 /// Wraps `val` into the range `[min, max]` using modular arithmetics.
 ///
 /// The range must not be empty.
+#[must_use]
 #[inline]
 pub fn wrap<T>(mut val: T, min: T, max: T) -> T
 where
@@ -219,6 +220,7 @@ where
 ///     * If `min < val < max`, this returns `val`.
 ///     * If `val <= min`, this returns `min`.
 ///     * If `val >= max`, this returns `max`.
+#[must_use]
 #[inline]
 pub fn clamp<T: PartialOrd>(val: T, min: T, max: T) -> T {
     if val > min {
diff --git a/src/linalg/bidiagonal.rs b/src/linalg/bidiagonal.rs
index a125ce18..6a462988 100644
--- a/src/linalg/bidiagonal.rs
+++ b/src/linalg/bidiagonal.rs
@@ -153,6 +153,7 @@ where
 
     /// Indicates whether this decomposition contains an upper-diagonal matrix.
     #[inline]
+    #[must_use]
     pub fn is_upper_diagonal(&self) -> bool {
         self.upper_diagonal
     }
@@ -188,6 +189,7 @@ where
 
     /// Retrieves the upper trapezoidal submatrix `R` of this decomposition.
     #[inline]
+    #[must_use]
     pub fn d(&self) -> OMatrix<T, DimMinimum<R, C>, DimMinimum<R, C>>
     where
         DefaultAllocator: Allocator<T, DimMinimum<R, C>, DimMinimum<R, C>>,
@@ -207,6 +209,7 @@ where
     /// Computes the orthogonal matrix `U` of this `U * D * V` decomposition.
     // TODO: code duplication with householder::assemble_q.
     // Except that we are returning a rectangular matrix here.
+    #[must_use]
     pub fn u(&self) -> OMatrix<T, R, DimMinimum<R, C>>
     where
         DefaultAllocator: Allocator<T, R, DimMinimum<R, C>>,
@@ -237,6 +240,7 @@ where
     }
 
     /// Computes the orthogonal matrix `V_t` of this `U * D * V_t` decomposition.
+    #[must_use]
     pub fn v_t(&self) -> OMatrix<T, DimMinimum<R, C>, C>
     where
         DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>,
@@ -274,6 +278,7 @@ where
     }
 
     /// The diagonal part of this decomposed matrix.
+    #[must_use]
     pub fn diagonal(&self) -> OVector<T::RealField, DimMinimum<R, C>>
     where
         DefaultAllocator: Allocator<T::RealField, DimMinimum<R, C>>,
@@ -282,6 +287,7 @@ where
     }
 
     /// The off-diagonal part of this decomposed matrix.
+    #[must_use]
     pub fn off_diagonal(&self) -> OVector<T::RealField, DimDiff<DimMinimum<R, C>, U1>>
     where
         DefaultAllocator: Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>,
diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs
index 24a6bcf5..f66fb42f 100644
--- a/src/linalg/cholesky.rs
+++ b/src/linalg/cholesky.rs
@@ -92,6 +92,7 @@ where
 
     /// Retrieves the lower-triangular factor of the Cholesky decomposition with its strictly
     /// uppen-triangular part filled with zeros.
+    #[must_use]
     pub fn l(&self) -> OMatrix<T, D, D> {
         self.chol.lower_triangle()
     }
@@ -101,6 +102,7 @@ where
     ///
     /// This is an allocation-less version of `self.l()`. The values of the strict upper-triangular
     /// part are garbage and should be ignored by further computations.
+    #[must_use]
     pub fn l_dirty(&self) -> &OMatrix<T, D, D> {
         &self.chol
     }
@@ -119,6 +121,7 @@ where
 
     /// Returns the solution of the system `self * x = b` where `self` is the decomposed matrix and
     /// `x` the unknown.
+    #[must_use = "Did you mean to use solve_mut()?"]
     pub fn solve<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<T, R2, C2, S2>) -> OMatrix<T, R2, C2>
     where
         S2: Storage<T, R2, C2>,
@@ -131,6 +134,7 @@ where
     }
 
     /// Computes the inverse of the decomposed matrix.
+    #[must_use]
     pub fn inverse(&self) -> OMatrix<T, D, D> {
         let shape = self.chol.data.shape();
         let mut res = OMatrix::identity_generic(shape.0, shape.1);
@@ -140,6 +144,7 @@ where
     }
 
     /// Computes the determinant of the decomposed matrix.
+    #[must_use]
     pub fn determinant(&self) -> T::SimdRealField {
         let dim = self.chol.nrows();
         let mut prod_diag = T::one();
@@ -287,6 +292,7 @@ where
 
     /// Updates the decomposition such that we get the decomposition of the factored matrix with its `j`th column removed.
     /// Since the matrix is square, the `j`th row will also be removed.
+    #[must_use]
     pub fn remove_column(&self, j: usize) -> Cholesky<T, DimDiff<D, U1>>
     where
         D: DimSub<U1>,
diff --git a/src/linalg/col_piv_qr.rs b/src/linalg/col_piv_qr.rs
index 10997b43..1a56d2cb 100644
--- a/src/linalg/col_piv_qr.rs
+++ b/src/linalg/col_piv_qr.rs
@@ -95,6 +95,7 @@ where
 
     /// Retrieves the upper trapezoidal submatrix `R` of this decomposition.
     #[inline]
+    #[must_use]
     pub fn r(&self) -> OMatrix<T, DimMinimum<R, C>, C>
     where
         DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>,
@@ -126,6 +127,7 @@ where
     }
 
     /// Computes the orthogonal matrix `Q` of this decomposition.
+    #[must_use]
     pub fn q(&self) -> OMatrix<T, R, DimMinimum<R, C>>
     where
         DefaultAllocator: Allocator<T, R, DimMinimum<R, C>>,
@@ -150,6 +152,7 @@ where
     }
     /// Retrieves the column permutation of this decomposition.
     #[inline]
+    #[must_use]
     pub fn p(&self) -> &PermutationSequence<DimMinimum<R, C>> {
         &self.p
     }
@@ -201,6 +204,7 @@ where
     /// Solves the linear system `self * x = b`, where `x` is the unknown to be determined.
     ///
     /// Returns `None` if `self` is not invertible.
+    #[must_use = "Did you mean to use solve_mut()?"]
     pub fn solve<R2: Dim, C2: Dim, S2>(
         &self,
         b: &Matrix<T, R2, C2, S2>,
@@ -283,6 +287,7 @@ where
     /// Computes the inverse of the decomposed matrix.
     ///
     /// Returns `None` if the decomposed matrix is not invertible.
+    #[must_use]
     pub fn try_inverse(&self) -> Option<OMatrix<T, D, D>> {
         assert!(
             self.col_piv_qr.is_square(),
@@ -301,6 +306,7 @@ where
     }
 
     /// Indicates if the decomposed matrix is invertible.
+    #[must_use]
     pub fn is_invertible(&self) -> bool {
         assert!(
             self.col_piv_qr.is_square(),
@@ -317,6 +323,7 @@ where
     }
 
     /// Computes the determinant of the decomposed matrix.
+    #[must_use]
     pub fn determinant(&self) -> T {
         let dim = self.col_piv_qr.nrows();
         assert!(
diff --git a/src/linalg/convolution.rs b/src/linalg/convolution.rs
index 2729b66b..36cea3a0 100644
--- a/src/linalg/convolution.rs
+++ b/src/linalg/convolution.rs
@@ -112,6 +112,7 @@ impl<T: RealField, D1: Dim, S1: Storage<T, D1>> Vector<T, D1, S1> {
     ///
     /// # Errors
     /// Inputs must satisfy `self.len() >= kernel.len() > 0`.
+    #[must_use]
     pub fn convolve_same<D2, S2>(&self, kernel: Vector<T, D2, S2>) -> OVector<T, D1>
     where
         D2: Dim,
diff --git a/src/linalg/determinant.rs b/src/linalg/determinant.rs
index aa04ff3f..22b681f5 100644
--- a/src/linalg/determinant.rs
+++ b/src/linalg/determinant.rs
@@ -12,6 +12,7 @@ impl<T: ComplexField, D: DimMin<D, Output = D>, S: Storage<T, D, D>> SquareMatri
     ///
     /// If the matrix has a dimension larger than 3, an LU decomposition is used.
     #[inline]
+    #[must_use]
     pub fn determinant(&self) -> T
     where
         DefaultAllocator: Allocator<T, D, D> + Allocator<(usize, usize), D>,
diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index 598c6350..acb2b9e3 100644
--- a/src/linalg/exp.rs
+++ b/src/linalg/exp.rs
@@ -435,6 +435,7 @@ where
         + Allocator<T::RealField, D, D>,
 {
     /// Computes exponential of this matrix
+    #[must_use]
     pub fn exp(&self) -> Self {
         // Simple case
         if self.nrows() == 1 {
diff --git a/src/linalg/full_piv_lu.rs b/src/linalg/full_piv_lu.rs
index 1555913b..f08af55c 100644
--- a/src/linalg/full_piv_lu.rs
+++ b/src/linalg/full_piv_lu.rs
@@ -96,6 +96,7 @@ where
 
     /// The lower triangular matrix of this decomposition.
     #[inline]
+    #[must_use]
     pub fn l(&self) -> OMatrix<T, R, DimMinimum<R, C>>
     where
         DefaultAllocator: Allocator<T, R, DimMinimum<R, C>>,
@@ -109,6 +110,7 @@ where
 
     /// The upper triangular matrix of this decomposition.
     #[inline]
+    #[must_use]
     pub fn u(&self) -> OMatrix<T, DimMinimum<R, C>, C>
     where
         DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>,
@@ -119,12 +121,14 @@ where
 
     /// The row permutations of this decomposition.
     #[inline]
+    #[must_use]
     pub fn p(&self) -> &PermutationSequence<DimMinimum<R, C>> {
         &self.p
     }
 
     /// The column permutations of this decomposition.
     #[inline]
+    #[must_use]
     pub fn q(&self) -> &PermutationSequence<DimMinimum<R, C>> {
         &self.q
     }
@@ -159,6 +163,7 @@ where
     /// Solves the linear system `self * x = b`, where `x` is the unknown to be determined.
     ///
     /// Returns `None` if the decomposed matrix is not invertible.
+    #[must_use = "Did you mean to use solve_mut()?"]
     pub fn solve<R2: Dim, C2: Dim, S2>(
         &self,
         b: &Matrix<T, R2, C2, S2>,
@@ -210,6 +215,7 @@ where
     /// Computes the inverse of the decomposed matrix.
     ///
     /// Returns `None` if the decomposed matrix is not invertible.
+    #[must_use]
     pub fn try_inverse(&self) -> Option<OMatrix<T, D, D>> {
         assert!(
             self.lu.is_square(),
@@ -227,6 +233,7 @@ where
     }
 
     /// Indicates if the decomposed matrix is invertible.
+    #[must_use]
     pub fn is_invertible(&self) -> bool {
         assert!(
             self.lu.is_square(),
@@ -238,6 +245,7 @@ where
     }
 
     /// Computes the determinant of the decomposed matrix.
+    #[must_use]
     pub fn determinant(&self) -> T {
         assert!(
             self.lu.is_square(),
diff --git a/src/linalg/givens.rs b/src/linalg/givens.rs
index 6073634e..8be91fe1 100644
--- a/src/linalg/givens.rs
+++ b/src/linalg/givens.rs
@@ -89,16 +89,19 @@ impl<T: ComplexField> GivensRotation<T> {
     }
 
     /// The cos part of this roration.
+    #[must_use]
     pub fn c(&self) -> T::RealField {
         self.c
     }
 
     /// The sin part of this roration.
+    #[must_use]
     pub fn s(&self) -> T {
         self.s
     }
 
     /// The inverse of this givens rotation.
+    #[must_use = "This function does not mutate self."]
     pub fn inverse(&self) -> Self {
         Self {
             c: self.c,
diff --git a/src/linalg/hessenberg.rs b/src/linalg/hessenberg.rs
index be6f667c..6b8ecfee 100644
--- a/src/linalg/hessenberg.rs
+++ b/src/linalg/hessenberg.rs
@@ -116,6 +116,7 @@ where
     ///
     /// This is less efficient than `.unpack_h()` as it allocates a new matrix.
     #[inline]
+    #[must_use]
     pub fn h(&self) -> OMatrix<T, D, D> {
         let dim = self.hess.nrows();
         let mut res = self.hess.clone();
@@ -126,6 +127,7 @@ where
     }
 
     /// Computes the orthogonal matrix `Q` of this decomposition.
+    #[must_use]
     pub fn q(&self) -> OMatrix<T, D, D> {
         householder::assemble_q(&self.hess, self.subdiag.as_slice())
     }
diff --git a/src/linalg/lu.rs b/src/linalg/lu.rs
index fce3f3cb..36a00807 100644
--- a/src/linalg/lu.rs
+++ b/src/linalg/lu.rs
@@ -127,6 +127,7 @@ where
 
     /// The lower triangular matrix of this decomposition.
     #[inline]
+    #[must_use]
     pub fn l(&self) -> OMatrix<T, R, DimMinimum<R, C>>
     where
         DefaultAllocator: Allocator<T, R, DimMinimum<R, C>>,
@@ -170,6 +171,7 @@ where
 
     /// The upper triangular matrix of this decomposition.
     #[inline]
+    #[must_use]
     pub fn u(&self) -> OMatrix<T, DimMinimum<R, C>, C>
     where
         DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>,
@@ -180,6 +182,7 @@ where
 
     /// The row permutations of this decomposition.
     #[inline]
+    #[must_use]
     pub fn p(&self) -> &PermutationSequence<DimMinimum<R, C>> {
         &self.p
     }
@@ -213,6 +216,7 @@ where
     /// Solves the linear system `self * x = b`, where `x` is the unknown to be determined.
     ///
     /// Returns `None` if `self` is not invertible.
+    #[must_use = "Did you mean to use solve_mut()?"]
     pub fn solve<R2: Dim, C2: Dim, S2>(
         &self,
         b: &Matrix<T, R2, C2, S2>,
@@ -257,6 +261,7 @@ where
     /// Computes the inverse of the decomposed matrix.
     ///
     /// Returns `None` if the matrix is not invertible.
+    #[must_use]
     pub fn try_inverse(&self) -> Option<OMatrix<T, D, D>> {
         assert!(
             self.lu.is_square(),
@@ -291,6 +296,7 @@ where
     }
 
     /// Computes the determinant of the decomposed matrix.
+    #[must_use]
     pub fn determinant(&self) -> T {
         let dim = self.lu.nrows();
         assert!(
@@ -307,6 +313,7 @@ where
     }
 
     /// Indicates if the decomposed matrix is invertible.
+    #[must_use]
     pub fn is_invertible(&self) -> bool {
         assert!(
             self.lu.is_square(),
diff --git a/src/linalg/permutation_sequence.rs b/src/linalg/permutation_sequence.rs
index 4cc8cb25..ea868b5a 100644
--- a/src/linalg/permutation_sequence.rs
+++ b/src/linalg/permutation_sequence.rs
@@ -140,17 +140,20 @@ where
     }
 
     /// The number of non-identity permutations applied by this sequence.
+    #[must_use]
     pub fn len(&self) -> usize {
         self.len
     }
 
     /// Returns true if the permutation sequence contains no elements.
+    #[must_use]
     pub fn is_empty(&self) -> bool {
         self.len() == 0
     }
 
     /// The determinant of the matrix corresponding to this permutation.
     #[inline]
+    #[must_use]
     pub fn determinant<T: One + ClosedNeg>(&self) -> T {
         if self.len % 2 == 0 {
             T::one()
diff --git a/src/linalg/qr.rs b/src/linalg/qr.rs
index 6d5274e4..4bdbb364 100644
--- a/src/linalg/qr.rs
+++ b/src/linalg/qr.rs
@@ -70,6 +70,7 @@ where
 
     /// Retrieves the upper trapezoidal submatrix `R` of this decomposition.
     #[inline]
+    #[must_use]
     pub fn r(&self) -> OMatrix<T, DimMinimum<R, C>, C>
     where
         DefaultAllocator: Allocator<T, DimMinimum<R, C>, C>,
@@ -96,6 +97,7 @@ where
     }
 
     /// Computes the orthogonal matrix `Q` of this decomposition.
+    #[must_use]
     pub fn q(&self) -> OMatrix<T, R, DimMinimum<R, C>>
     where
         DefaultAllocator: Allocator<T, R, DimMinimum<R, C>>,
@@ -164,6 +166,7 @@ where
     /// Solves the linear system `self * x = b`, where `x` is the unknown to be determined.
     ///
     /// Returns `None` if `self` is not invertible.
+    #[must_use = "Did you mean to use solve_mut()?"]
     pub fn solve<R2: Dim, C2: Dim, S2>(
         &self,
         b: &Matrix<T, R2, C2, S2>,
@@ -243,6 +246,7 @@ where
     /// Computes the inverse of the decomposed matrix.
     ///
     /// Returns `None` if the decomposed matrix is not invertible.
+    #[must_use]
     pub fn try_inverse(&self) -> Option<OMatrix<T, D, D>> {
         assert!(
             self.qr.is_square(),
@@ -261,6 +265,7 @@ where
     }
 
     /// Indicates if the decomposed matrix is invertible.
+    #[must_use]
     pub fn is_invertible(&self) -> bool {
         assert!(
             self.qr.is_square(),
diff --git a/src/linalg/schur.rs b/src/linalg/schur.rs
index 1903b717..16b8c66d 100644
--- a/src/linalg/schur.rs
+++ b/src/linalg/schur.rs
@@ -385,6 +385,7 @@ where
     /// Computes the real eigenvalues of the decomposed matrix.
     ///
     /// Return `None` if some eigenvalues are complex.
+    #[must_use]
     pub fn eigenvalues(&self) -> Option<OVector<T, D>> {
         let mut out = unsafe {
             crate::unimplemented_or_uninitialized_generic!(self.t.data.shape().0, Const::<1>)
@@ -397,6 +398,7 @@ where
     }
 
     /// Computes the complex eigenvalues of the decomposed matrix.
+    #[must_use]
     pub fn complex_eigenvalues(&self) -> OVector<NumComplex<T>, D>
     where
         T: RealField,
@@ -509,6 +511,7 @@ where
         + Allocator<T, D>,
 {
     /// Computes the eigenvalues of this matrix.
+    #[must_use]
     pub fn eigenvalues(&self) -> Option<OVector<T, D>> {
         assert!(
             self.is_square(),
@@ -551,6 +554,7 @@ where
     }
 
     /// Computes the eigenvalues of this matrix.
+    #[must_use]
     pub fn complex_eigenvalues(&self) -> OVector<NumComplex<T>, D>
     // TODO: add balancing?
     where
diff --git a/src/linalg/solve.rs b/src/linalg/solve.rs
index 8fdab47e..7f9b7dae 100644
--- a/src/linalg/solve.rs
+++ b/src/linalg/solve.rs
@@ -10,6 +10,7 @@ use crate::base::{DVectorSlice, DefaultAllocator, Matrix, OMatrix, SquareMatrix,
 impl<T: ComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
     /// Computes the solution of the linear system `self . x = b` where `x` is the unknown and only
     /// the lower-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[must_use = "Did you mean to use solve_lower_triangular_mut()?"]
     #[inline]
     pub fn solve_lower_triangular<R2: Dim, C2: Dim, S2>(
         &self,
@@ -30,6 +31,7 @@ impl<T: ComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
 
     /// Computes the solution of the linear system `self . x = b` where `x` is the unknown and only
     /// the upper-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[must_use = "Did you mean to use solve_upper_triangular_mut()?"]
     #[inline]
     pub fn solve_upper_triangular<R2: Dim, C2: Dim, S2>(
         &self,
@@ -186,6 +188,7 @@ impl<T: ComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
      */
     /// Computes the solution of the linear system `self.transpose() . x = b` where `x` is the unknown and only
     /// the lower-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[must_use = "Did you mean to use tr_solve_lower_triangular_mut()?"]
     #[inline]
     pub fn tr_solve_lower_triangular<R2: Dim, C2: Dim, S2>(
         &self,
@@ -206,6 +209,7 @@ impl<T: ComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
 
     /// Computes the solution of the linear system `self.transpose() . x = b` where `x` is the unknown and only
     /// the upper-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[must_use = "Did you mean to use tr_solve_upper_triangular_mut()?"]
     #[inline]
     pub fn tr_solve_upper_triangular<R2: Dim, C2: Dim, S2>(
         &self,
@@ -276,6 +280,7 @@ impl<T: ComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
 
     /// Computes the solution of the linear system `self.adjoint() . x = b` where `x` is the unknown and only
     /// the lower-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[must_use = "Did you mean to use ad_solve_lower_triangular_mut()?"]
     #[inline]
     pub fn ad_solve_lower_triangular<R2: Dim, C2: Dim, S2>(
         &self,
@@ -296,6 +301,7 @@ impl<T: ComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
 
     /// Computes the solution of the linear system `self.adjoint() . x = b` where `x` is the unknown and only
     /// the upper-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[must_use = "Did you mean to use ad_solve_upper_triangular_mut()?"]
     #[inline]
     pub fn ad_solve_upper_triangular<R2: Dim, C2: Dim, S2>(
         &self,
@@ -443,6 +449,7 @@ impl<T: ComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
 impl<T: SimdComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
     /// Computes the solution of the linear system `self . x = b` where `x` is the unknown and only
     /// the lower-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[must_use = "Did you mean to use solve_lower_triangular_unchecked_mut()?"]
     #[inline]
     pub fn solve_lower_triangular_unchecked<R2: Dim, C2: Dim, S2>(
         &self,
@@ -460,6 +467,7 @@ impl<T: SimdComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
 
     /// Computes the solution of the linear system `self . x = b` where `x` is the unknown and only
     /// the upper-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[must_use = "Did you mean to use solve_upper_triangular_unchecked_mut()?"]
     #[inline]
     pub fn solve_upper_triangular_unchecked<R2: Dim, C2: Dim, S2>(
         &self,
@@ -578,6 +586,7 @@ impl<T: SimdComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
      */
     /// Computes the solution of the linear system `self.transpose() . x = b` where `x` is the unknown and only
     /// the lower-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[must_use = "Did you mean to use tr_solve_lower_triangular_unchecked_mut()?"]
     #[inline]
     pub fn tr_solve_lower_triangular_unchecked<R2: Dim, C2: Dim, S2>(
         &self,
@@ -595,6 +604,7 @@ impl<T: SimdComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
 
     /// Computes the solution of the linear system `self.transpose() . x = b` where `x` is the unknown and only
     /// the upper-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[must_use = "Did you mean to use tr_solve_upper_triangular_unchecked_mut()?"]
     #[inline]
     pub fn tr_solve_upper_triangular_unchecked<R2: Dim, C2: Dim, S2>(
         &self,
@@ -648,6 +658,7 @@ impl<T: SimdComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
 
     /// Computes the solution of the linear system `self.adjoint() . x = b` where `x` is the unknown and only
     /// the lower-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[must_use = "Did you mean to use ad_solve_lower_triangular_unchecked_mut()?"]
     #[inline]
     pub fn ad_solve_lower_triangular_unchecked<R2: Dim, C2: Dim, S2>(
         &self,
@@ -665,6 +676,7 @@ impl<T: SimdComplexField, D: Dim, S: Storage<T, D, D>> SquareMatrix<T, D, S> {
 
     /// Computes the solution of the linear system `self.adjoint() . x = b` where `x` is the unknown and only
     /// the upper-triangular part of `self` (including the diagonal) is considered not-zero.
+    #[must_use = "Did you mean to use ad_solve_upper_triangular_unchecked_mut()?"]
     #[inline]
     pub fn ad_solve_upper_triangular_unchecked<R2: Dim, C2: Dim, S2>(
         &self,
diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs
index 2f5e6a63..f559c9fa 100644
--- a/src/linalg/svd.rs
+++ b/src/linalg/svd.rs
@@ -502,6 +502,7 @@ where
 
     /// Computes the rank of the decomposed matrix, i.e., the number of singular values greater
     /// than `eps`.
+    #[must_use]
     pub fn rank(&self, eps: T::RealField) -> usize {
         assert!(
             eps >= T::RealField::zero(),
@@ -615,6 +616,7 @@ where
         + Allocator<T::RealField, DimDiff<DimMinimum<R, C>, U1>>,
 {
     /// Computes the singular values of this matrix.
+    #[must_use]
     pub fn singular_values(&self) -> OVector<T::RealField, DimMinimum<R, C>> {
         SVD::new(self.clone_owned(), false, false).singular_values
     }
@@ -622,6 +624,7 @@ where
     /// Computes the rank of this matrix.
     ///
     /// All singular values below `eps` are considered equal to 0.
+    #[must_use]
     pub fn rank(&self, eps: T::RealField) -> usize {
         let svd = SVD::new(self.clone_owned(), false, false);
         svd.rank(eps)
diff --git a/src/linalg/symmetric_eigen.rs b/src/linalg/symmetric_eigen.rs
index 6e3b67fe..25dab08c 100644
--- a/src/linalg/symmetric_eigen.rs
+++ b/src/linalg/symmetric_eigen.rs
@@ -268,6 +268,7 @@ where
     /// Rebuild the original matrix.
     ///
     /// This is useful if some of the eigenvalues have been manually modified.
+    #[must_use]
     pub fn recompose(&self) -> OMatrix<T, D, D> {
         let mut u_t = self.eigenvectors.clone();
         for i in 0..self.eigenvalues.len() {
@@ -311,6 +312,7 @@ where
     /// Computes the eigenvalues of this symmetric matrix.
     ///
     /// Only the lower-triangular part of the matrix is read.
+    #[must_use]
     pub fn symmetric_eigenvalues(&self) -> OVector<T::RealField, D> {
         SymmetricEigen::do_decompose(
             self.clone_owned(),
diff --git a/src/linalg/symmetric_tridiagonal.rs b/src/linalg/symmetric_tridiagonal.rs
index 9937f76b..c7e87ba8 100644
--- a/src/linalg/symmetric_tridiagonal.rs
+++ b/src/linalg/symmetric_tridiagonal.rs
@@ -131,6 +131,7 @@ where
     }
 
     /// The diagonal components of this decomposition.
+    #[must_use]
     pub fn diagonal(&self) -> OVector<T::RealField, D>
     where
         DefaultAllocator: Allocator<T::RealField, D>,
@@ -139,6 +140,7 @@ where
     }
 
     /// The off-diagonal components of this decomposition.
+    #[must_use]
     pub fn off_diagonal(&self) -> OVector<T::RealField, DimDiff<D, U1>>
     where
         DefaultAllocator: Allocator<T::RealField, DimDiff<D, U1>>,
@@ -147,6 +149,7 @@ where
     }
 
     /// Computes the orthogonal matrix `Q` of this decomposition.
+    #[must_use]
     pub fn q(&self) -> OMatrix<T, D, D> {
         householder::assemble_q(&self.tri, self.off_diagonal.as_slice())
     }
diff --git a/src/linalg/udu.rs b/src/linalg/udu.rs
index 4d3aca6c..7b4a9cc9 100644
--- a/src/linalg/udu.rs
+++ b/src/linalg/udu.rs
@@ -92,6 +92,7 @@ where
     }
 
     /// Returns the diagonal elements as a matrix
+    #[must_use]
     pub fn d_matrix(&self) -> OMatrix<T, D, D> {
         OMatrix::from_diagonal(&self.d)
     }
diff --git a/src/sparse/cs_matrix.rs b/src/sparse/cs_matrix.rs
index c4a0c901..74e13719 100644
--- a/src/sparse/cs_matrix.rs
+++ b/src/sparse/cs_matrix.rs
@@ -119,16 +119,19 @@ where
     DefaultAllocator: Allocator<usize, C>,
 {
     /// The value buffer of this storage.
+    #[must_use]
     pub fn values(&self) -> &[T] {
         &self.vals
     }
 
     /// The column shifts buffer.
+    #[must_use]
     pub fn p(&self) -> &[usize] {
         self.p.as_slice()
     }
 
     /// The row index buffers.
+    #[must_use]
     pub fn i(&self) -> &[usize] {
         &self.i
     }
@@ -356,27 +359,32 @@ impl<T: Scalar, R: Dim, C: Dim, S: CsStorage<T, R, C>> CsMatrix<T, R, C, S> {
     }
 
     /// The size of the data buffer.
+    #[must_use]
     pub fn len(&self) -> usize {
         self.data.len()
     }
 
     /// The number of rows of this matrix.
+    #[must_use]
     pub fn nrows(&self) -> usize {
         self.data.shape().0.value()
     }
 
     /// The number of rows of this matrix.
+    #[must_use]
     pub fn ncols(&self) -> usize {
         self.data.shape().1.value()
     }
 
     /// The shape of this matrix.
+    #[must_use]
     pub fn shape(&self) -> (usize, usize) {
         let (nrows, ncols) = self.data.shape();
         (nrows.value(), ncols.value())
     }
 
     /// Whether this matrix is square or not.
+    #[must_use]
     pub fn is_square(&self) -> bool {
         let (nrows, ncols) = self.data.shape();
         nrows.value() == ncols.value()
@@ -391,6 +399,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: CsStorage<T, R, C>> CsMatrix<T, R, C, S> {
     /// If at any time this `is_sorted` method returns `false`, then, something went wrong
     /// and an issue should be open on the nalgebra repository with details on how to reproduce
     /// this.
+    #[must_use]
     pub fn is_sorted(&self) -> bool {
         for j in 0..self.ncols() {
             let mut curr = None;
@@ -409,6 +418,7 @@ impl<T: Scalar, R: Dim, C: Dim, S: CsStorage<T, R, C>> CsMatrix<T, R, C, S> {
     }
 
     /// Computes the transpose of this sparse matrix.
+    #[must_use = "This function does not mutate the matrix. Consider using the return value or removing the function call. There's also transpose_mut() for square matrices."]
     pub fn transpose(&self) -> CsMatrix<T, C, R>
     where
         DefaultAllocator: Allocator<usize, R>,
diff --git a/src/sparse/cs_matrix_cholesky.rs b/src/sparse/cs_matrix_cholesky.rs
index 4cd61d2a..6d52d0a6 100644
--- a/src/sparse/cs_matrix_cholesky.rs
+++ b/src/sparse/cs_matrix_cholesky.rs
@@ -67,6 +67,7 @@ where
     }
 
     /// The lower-triangular matrix of the cholesky decomposition.
+    #[must_use]
     pub fn l(&self) -> Option<&CsMatrix<T, D, D>> {
         if self.ok {
             Some(&self.l)
diff --git a/src/sparse/cs_matrix_solve.rs b/src/sparse/cs_matrix_solve.rs
index 52065747..43c5c2c7 100644
--- a/src/sparse/cs_matrix_solve.rs
+++ b/src/sparse/cs_matrix_solve.rs
@@ -6,6 +6,7 @@ use crate::{Const, DefaultAllocator, Dim, Matrix, OMatrix, OVector, RealField};
 
 impl<T: RealField, D: Dim, S: CsStorage<T, D, D>> CsMatrix<T, D, D, S> {
     /// Solve a lower-triangular system with a dense right-hand-side.
+    #[must_use = "Did you mean to use solve_lower_triangular_mut()?"]
     pub fn solve_lower_triangular<R2: Dim, C2: Dim, S2>(
         &self,
         b: &Matrix<T, R2, C2, S2>,
@@ -24,6 +25,7 @@ impl<T: RealField, D: Dim, S: CsStorage<T, D, D>> CsMatrix<T, D, D, S> {
     }
 
     /// Solve a lower-triangular system with `self` transposed and a dense right-hand-side.
+    #[must_use = "Did you mean to use tr_solve_lower_triangular_mut()?"]
     pub fn tr_solve_lower_triangular<R2: Dim, C2: Dim, S2>(
         &self,
         b: &Matrix<T, R2, C2, S2>,
@@ -135,6 +137,7 @@ impl<T: RealField, D: Dim, S: CsStorage<T, D, D>> CsMatrix<T, D, D, S> {
     }
 
     /// Solve a lower-triangular system with a sparse right-hand-side.
+    #[must_use]
     pub fn solve_lower_triangular_cs<D2: Dim, S2>(
         &self,
         b: &CsVector<T, D2, S2>,