diff --git a/nalgebra-lapack/src/cholesky.rs b/nalgebra-lapack/src/cholesky.rs
index 4b3c8a39..1bddbc20 100644
--- a/nalgebra-lapack/src/cholesky.rs
+++ b/nalgebra-lapack/src/cholesky.rs
@@ -1,8 +1,9 @@
 use num::Zero;
 use num_complex::Complex;
 
-use na::{Scalar, DefaultAllocator, MatrixN, MatrixMN};
+use na::{Scalar, DefaultAllocator, Matrix, MatrixN, MatrixMN};
 use na::dimension::Dim;
+use na::storage::Storage;
 use na::allocator::Allocator;
 
 use lapack::fortran as interface;
@@ -48,8 +49,22 @@ impl<N: CholeskyScalar + Zero, D: Dim> Cholesky<N, D>
 
     /// Solves the symmetric-definite-positive linear system `self * x = b`, where `x` is the
     /// unknown to be determined.
-    pub fn solve<R2: Dim, C2: Dim>(&self, mut b: MatrixMN<N, R2, C2>)
-        -> Option<MatrixMN<N, R2, C2>>
+    pub fn solve<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<N, R2, C2, S2>) -> Option<MatrixMN<N, R2, C2>>
+        where S2: Storage<N, R2, C2>,
+              DefaultAllocator: Allocator<N, R2, C2> {
+
+        let mut res = b.clone_owned();
+        if self.solve_mut(&mut res) {
+            Some(res)
+        }
+        else {
+            None
+        }
+    }
+
+    /// Solves in-place the symmetric-definite-positive linear system `self * x = b`, where `x` is
+    /// the unknown to be determined.
+    pub fn solve_mut<R2: Dim, C2: Dim>(&self, b: &mut MatrixMN<N, R2, C2>) -> bool
         where DefaultAllocator: Allocator<N, R2, C2> {
 
         let dim = self.l.nrows();
@@ -62,9 +77,7 @@ impl<N: CholeskyScalar + Zero, D: Dim> Cholesky<N, D>
         let mut info = 0;
 
         N::xpotrs(b'L', dim as i32, nrhs, self.l.as_slice(), lda, b.as_mut_slice(), ldb, &mut info);
-        lapack_check!(info);
-
-        Some(b)
+        lapack_test!(info)
     }
 
     /// Computes the inverse of the decomposed matrix.
diff --git a/nalgebra-lapack/src/lapack_check.rs b/nalgebra-lapack/src/lapack_check.rs
index ee777706..217699dc 100644
--- a/nalgebra-lapack/src/lapack_check.rs
+++ b/nalgebra-lapack/src/lapack_check.rs
@@ -19,3 +19,9 @@ macro_rules! lapack_panic(
         assert!($info == 0, "Lapack error.");
     );
 );
+
+macro_rules! lapack_test(
+    ($info: expr) => (
+        $info == 0
+    );
+);
diff --git a/nalgebra-lapack/src/lu.rs b/nalgebra-lapack/src/lu.rs
index cc3337b0..b7655715 100644
--- a/nalgebra-lapack/src/lu.rs
+++ b/nalgebra-lapack/src/lu.rs
@@ -105,8 +105,7 @@ impl<N: LUScalar, R: Dim, C: Dim> LU<N, R, C>
                   1, self.p.len() as i32, self.p.as_slice(), -1);
     }
 
-    fn generic_solve<R2: Dim, C2: Dim>(&self, trans: u8, mut b: MatrixMN<N, R2, C2>)
-        -> Option<MatrixMN<N, R2, C2>>
+    fn generic_solve_mut<R2: Dim, C2: Dim>(&self, trans: u8, b: &mut MatrixMN<N, R2, C2>) -> bool
         where DefaultAllocator: Allocator<N, R2, C2> +
                                 Allocator<i32, R2> {
 
@@ -121,39 +120,89 @@ impl<N: LUScalar, R: Dim, C: Dim> LU<N, R, C>
         let mut info = 0;
 
         N::xgetrs(trans, dim as i32, nrhs, self.lu.as_slice(), lda, self.p.as_slice(),
-                      b.as_mut_slice(), ldb, &mut info);
-        lapack_check!(info);
-
-        Some(b)
+                  b.as_mut_slice(), ldb, &mut info);
+        lapack_test!(info)
     }
 
     /// Solves the linear system `self * x = b`, where `x` is the unknown to be determined.
-    pub fn solve<R2: Dim, C2: Dim>(&self, b: MatrixMN<N, R2, C2>)
-        -> Option<MatrixMN<N, R2, C2>>
-        where DefaultAllocator: Allocator<N, R2, C2> +
+    pub fn solve<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<N, R2, C2, S2>) -> Option<MatrixMN<N, R2, C2>>
+        where S2: Storage<N, R2, C2>,
+              DefaultAllocator: Allocator<N, R2, C2> +
                                 Allocator<i32, R2> {
 
-        self.generic_solve(b'N', b)
+        let mut res = b.clone_owned();
+        if self.generic_solve_mut(b'N', &mut res) {
+            Some(res)
+        }
+        else {
+            None
+        }
     }
 
     /// Solves the linear system `self.transpose() * x = b`, where `x` is the unknown to be
     /// determined.
-    pub fn solve_transpose<R2: Dim, C2: Dim>(&self, b: MatrixMN<N, R2, C2>)
+    pub fn solve_transpose<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<N, R2, C2, S2>)
         -> Option<MatrixMN<N, R2, C2>>
-        where DefaultAllocator: Allocator<N, R2, C2> +
+        where S2: Storage<N, R2, C2>,
+              DefaultAllocator: Allocator<N, R2, C2> +
                                 Allocator<i32, R2> {
 
-        self.generic_solve(b'T', b)
+        let mut res = b.clone_owned();
+        if self.generic_solve_mut(b'T', &mut res) {
+            Some(res)
+        }
+        else {
+            None
+        }
     }
 
     /// Solves the linear system `self.conjugate_transpose() * x = b`, where `x` is the unknown to
     /// be determined.
-    pub fn solve_conjugate_transpose<R2: Dim, C2: Dim>(&self, b: MatrixMN<N, R2, C2>)
+    pub fn solve_conjugate_transpose<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<N, R2, C2, S2>)
         -> Option<MatrixMN<N, R2, C2>>
+        where S2: Storage<N, R2, C2>,
+              DefaultAllocator: Allocator<N, R2, C2> +
+                                Allocator<i32, R2> {
+
+        let mut res = b.clone_owned();
+        if self.generic_solve_mut(b'T', &mut res) {
+            Some(res)
+        }
+        else {
+            None
+        }
+    }
+
+    /// Solves in-place the linear system `self * x = b`, where `x` is the unknown to be determined.
+    ///
+    /// Retuns `false` if no solution was found (the decomposed matrix is singular).
+    pub fn solve_mut<R2: Dim, C2: Dim>(&self, b: &mut MatrixMN<N, R2, C2>) -> bool
         where DefaultAllocator: Allocator<N, R2, C2> +
                                 Allocator<i32, R2> {
 
-        self.generic_solve(b'T', b)
+        self.generic_solve_mut(b'N', b)
+    }
+
+    /// Solves in-place the linear system `self.transpose() * x = b`, where `x` is the unknown to be
+    /// determined.
+    ///
+    /// Retuns `false` if no solution was found (the decomposed matrix is singular).
+    pub fn solve_transpose_mut<R2: Dim, C2: Dim>(&self, b: &mut MatrixMN<N, R2, C2>) -> bool
+        where DefaultAllocator: Allocator<N, R2, C2> +
+                                Allocator<i32, R2> {
+
+        self.generic_solve_mut(b'T', b)
+    }
+
+    /// Solves in-place the linear system `self.conjugate_transpose() * x = b`, where `x` is the unknown to
+    /// be determined.
+    ///
+    /// Retuns `false` if no solution was found (the decomposed matrix is singular).
+    pub fn solve_conjugate_transpose_mut<R2: Dim, C2: Dim>(&self, b: &mut MatrixMN<N, R2, C2>) -> bool
+        where DefaultAllocator: Allocator<N, R2, C2> +
+                                Allocator<i32, R2> {
+
+        self.generic_solve_mut(b'T', b)
     }
 }
 
diff --git a/src/linalg/lu.rs b/src/linalg/lu.rs
index 770e64c7..532c6e32 100644
--- a/src/linalg/lu.rs
+++ b/src/linalg/lu.rs
@@ -169,7 +169,7 @@ impl<N: Real, D: DimMin<D, Output = D>> LU<N, D, D>
     ///
     /// Returns `None` if `self` is not invertible.
     pub fn solve<R2: Dim, C2: Dim, S2>(&self, b: &Matrix<N, R2, C2, S2>) -> Option<MatrixMN<N, R2, C2>>
-        where S2: StorageMut<N, R2, C2>,
+        where S2: Storage<N, R2, C2>,
               ShapeConstraint: SameNumberOfRows<R2, D>,
               DefaultAllocator: Allocator<N, R2, C2> {
         let mut res = b.clone_owned();