diff --git a/src/base/norm.rs b/src/base/norm.rs
index 251bf661..ab44ed02 100644
--- a/src/base/norm.rs
+++ b/src/base/norm.rs
@@ -340,6 +340,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`.
+    #[must_use = "This function does not mutate self but returns a new clamped version."]
     #[inline]
     pub fn cap_magnitude(&self, max: T::RealField) -> OMatrix<T, R, C>
     where
@@ -356,6 +357,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`.
+    #[must_use = "This function does not mutate self but returns a new clamped version."]
     #[inline]
     pub fn simd_cap_magnitude(&self, max: T::SimdRealField) -> OMatrix<T, R, C>
     where
diff --git a/src/lib.rs b/src/lib.rs
index 8378b272..f0f62fe1 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 = "This function does not mutate `val`."]
 #[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 = "This function does not mutate `val`."]
 #[inline]
 pub fn clamp<T: PartialOrd>(val: T, min: T, max: T) -> T {
     if val > min {
diff --git a/src/linalg/cholesky.rs b/src/linalg/cholesky.rs
index 24a6bcf5..eb2ebc9f 100644
--- a/src/linalg/cholesky.rs
+++ b/src/linalg/cholesky.rs
@@ -119,6 +119,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 +132,7 @@ where
     }
 
     /// Computes the inverse of the decomposed matrix.
+    #[must_use = "This function does not mutate self. Consider using the return value or dropping the function call."]
     pub fn inverse(&self) -> OMatrix<T, D, D> {
         let shape = self.chol.data.shape();
         let mut res = OMatrix::identity_generic(shape.0, shape.1);
diff --git a/src/linalg/col_piv_qr.rs b/src/linalg/col_piv_qr.rs
index 10997b43..3640785d 100644
--- a/src/linalg/col_piv_qr.rs
+++ b/src/linalg/col_piv_qr.rs
@@ -201,6 +201,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>,
diff --git a/src/linalg/exp.rs b/src/linalg/exp.rs
index 598c6350..b93055c5 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 = "This function does not mutate self. Consider using the return value or dropping the function call."]
     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..af9609de 100644
--- a/src/linalg/full_piv_lu.rs
+++ b/src/linalg/full_piv_lu.rs
@@ -159,6 +159,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>,
diff --git a/src/linalg/givens.rs b/src/linalg/givens.rs
index 6073634e..9a65f331 100644
--- a/src/linalg/givens.rs
+++ b/src/linalg/givens.rs
@@ -99,6 +99,7 @@ impl<T: ComplexField> GivensRotation<T> {
     }
 
     /// 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/lu.rs b/src/linalg/lu.rs
index fce3f3cb..855a53f7 100644
--- a/src/linalg/lu.rs
+++ b/src/linalg/lu.rs
@@ -213,6 +213,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>,
diff --git a/src/linalg/qr.rs b/src/linalg/qr.rs
index 6d5274e4..103916a0 100644
--- a/src/linalg/qr.rs
+++ b/src/linalg/qr.rs
@@ -164,6 +164,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>,
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/sparse/cs_matrix.rs b/src/sparse/cs_matrix.rs
index c4a0c901..2dc8467e 100644
--- a/src/sparse/cs_matrix.rs
+++ b/src/sparse/cs_matrix.rs
@@ -409,6 +409,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_solve.rs b/src/sparse/cs_matrix_solve.rs
index 52065747..33d811cf 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 = "This function has no side effects. Consider using the return value or removing the function call."]
     pub fn solve_lower_triangular_cs<D2: Dim, S2>(
         &self,
         b: &CsVector<T, D2, S2>,