diff --git a/collects/math/private/matrix/matrix-operations.rkt b/collects/math/private/matrix/matrix-operations.rkt
index b3ead31c02..24cfee1c2b 100644
--- a/collects/math/private/matrix/matrix-operations.rkt
+++ b/collects/math/private/matrix/matrix-operations.rkt
@@ -2,6 +2,7 @@
 
 (require racket/fixnum
          racket/list
+         racket/match
          math/array
          (only-in typed/racket conjugate)
          "../unsafe.rkt"
@@ -12,11 +13,12 @@
          "matrix-arithmetic.rkt"
          "matrix-basic.rkt"
          "matrix-column.rkt"
+         "utils.rkt"
          (for-syntax racket))
 
 ; TODO:
 ; 1. compute null space from QR factorization
-;    (better numerical stability than from Gauss elimnation)
+;    (better numerical stability than from Gauss elimination)
 ; 2. S+N decomposition
 ; 3. Linear least squares problems (data fitting)
 ; 4. Pseudo inverse
@@ -31,19 +33,16 @@
  matrix-swap-columns
  matrix-add-scaled-row
  ; reduction
- matrix-gauss-eliminate
- matrix-gauss-jordan-eliminate
- matrix-row-echelon-form
- matrix-reduced-row-echelon-form
+ matrix-gauss-elim
+ matrix-row-echelon
  ; invariant
  matrix-rank
  matrix-nullity
  matrix-determinant
- ; spaces
- ;matrix-column+null-space
+ matrix-determinant/row-reduction  ; for testing
+ matrix-invertible?
  ; solvers
  matrix-solve
- matrix-solve-many
  ; spaces
  matrix-column-space
  ; projection
@@ -98,7 +97,8 @@
            [(< j 0)
             (error 'matrix-scale-row "column index must be non-negative, got ~a" j)]
            [(not (< j (vector-ref ds 1)))
-            (error 'matrix-scale-row "column index must be smaller than the number of rows, got ~a" j)]
+            (error 'matrix-scale-row
+                   "column index must be smaller than the number of rows, got ~a" j)]
            [else
             (unsafe-build-array ds (λ: ([js : (Vectorof Index)]) 
                                      (if (= j (vector-ref js 1))
@@ -152,9 +152,11 @@
            [(< j 0)
             (error 'matrix-swap-columns "column index must be non-negative, got ~a" j)]
            [(not (< i (vector-ref ds 0)))
-            (error 'matrix-swap-columns "column index must be smaller than the number of columns, got ~a" i)]
+            (error 'matrix-swap-columns
+                   "column index must be smaller than the number of columns, got ~a" i)]
            [(not (< j (vector-ref ds 0)))
-            (error 'matrix-swap-columns "column index must be smaller than the number of columns, got ~a" j)]
+            (error 'matrix-swap-columns
+                   "column index must be smaller than the number of columns, got ~a" j)]
            [else
             (unsafe-build-array ds (λ: ([js : (Vectorof Index)]) 
                                      (cond
@@ -198,279 +200,226 @@
                                          (g js))))])))]))
 
 
-;;; GAUSS ELIMINATION / ROW ECHELON FORM
-
 (: unsafe-vector2d-ref (All (A) ((Vectorof (Vectorof A)) Index Index -> A)))
 (define (unsafe-vector2d-ref vss i j)
   (unsafe-vector-ref (unsafe-vector-ref vss i) j))
 
-(: find-partial-pivot (case-> ((Vectorof (Vectorof Real)) Index Index Index -> (U #f Index))
-                              ((Vectorof (Vectorof Number)) Index Index Index -> (U #f Index))))
+;; ===================================================================================================
+;; Gaussian elimination
+
+(: find-partial-pivot
+   (case-> ((Vectorof (Vectorof Real)) Index Index Index -> (Values Index Real))
+           ((Vectorof (Vectorof Number)) Index Index Index -> (Values Index Number))))
 ;; Find the element with maximum magnitude in a column
 (define (find-partial-pivot rows m i j)
-  (let loop ([#{l : Nonnegative-Fixnum} i] [#{max-current : Real} -inf.0] [#{max-index : Index} i])
+  (define l (fx+ i 1))
+  (define pivot (unsafe-vector2d-ref rows i j))
+  (define mag-pivot (magnitude pivot))
+  (let loop ([#{l : Nonnegative-Fixnum} l] [#{p : Index} i] [pivot pivot] [mag-pivot mag-pivot])
     (cond [(l . fx< . m)
-           (define v (magnitude (unsafe-vector2d-ref rows l j)))
-           (cond [(> v max-current)  (loop (fx+ l 1) v l)]
-                 [else  (loop (fx+ l 1) max-current max-index)])]
-          [else  max-index])))
+           (define new-pivot (unsafe-vector2d-ref rows l j))
+           (define mag-new-pivot (magnitude new-pivot))
+           (cond [(mag-new-pivot . > . mag-pivot)  (loop (fx+ l 1) l new-pivot mag-new-pivot)]
+                 [else  (loop (fx+ l 1) p pivot mag-pivot)])]
+          [else  (values p pivot)])))
 
-(: find-pivot (case-> ((Vectorof (Vectorof Real)) Index Index Index -> (U #f Index))
-                      ((Vectorof (Vectorof Number)) Index Index Index -> (U #f Index))))
-;; Find a non-zero element in a column
-(define (find-pivot rows m i j)
-  (let loop ([#{l : Nonnegative-Fixnum} i])
-    (cond [(l . fx>= . m)  #f]
-          [(not (zero? (unsafe-vector2d-ref rows l j)))  l]
-          [else  (loop (fx+ l 1))])))
+(: elim-rows!
+   (case-> ((Vectorof (Vectorof Real)) Index Index Index Real Nonnegative-Fixnum -> Void)
+           ((Vectorof (Vectorof Number)) Index Index Index Number Nonnegative-Fixnum -> Void)))
+(define (elim-rows! rows m i j pivot start)
+  (let loop ([#{l : Nonnegative-Fixnum} start])
+    (when (l . fx< . m)
+      (unless (l . fx= . i)
+        (define x_lj (unsafe-vector2d-ref rows l j))
+        (unless (zero? x_lj)
+          (vector-scaled-add! (unsafe-vector-ref rows l)
+                              (unsafe-vector-ref rows i)
+                              (- (/ x_lj pivot)))))
+      (loop (fx+ l 1)))))
 
-(: matrix-gauss-eliminate : 
-   (case-> ((Matrix Real)                 -> (Values (Matrix Real) (Listof Index)))
-           ((Matrix Real) Boolean         -> (Values (Matrix Real) (Listof Index)))
-           ((Matrix Real) Boolean Boolean -> (Values (Matrix Real) (Listof Index)))
-           ((Matrix Number)                 -> (Values (Matrix Number) (Listof Index)))
-           ((Matrix Number) Boolean         -> (Values (Matrix Number) (Listof Index)))
-           ((Matrix Number) Boolean Boolean -> (Values (Matrix Number) (Listof Index)))))
-;; Returns the result of Gaussian elimination and a list of column indexes that had no pivot value
-;; If `reduced?' is #t, the result is in *reduced* row-echelon form, and is unique (up to
-;; floating-point error)
-;; If `partial-pivoting?' is #t, the largest value in each column is used as the pivot
-(define (matrix-gauss-eliminate M [reduced? #f] [partial-pivoting? #t])
+(: matrix-gauss-elim (case-> ((Matrix Real)         -> (Values (Matrix Real) (Listof Index)))
+                             ((Matrix Real) Any     -> (Values (Matrix Real) (Listof Index)))
+                             ((Matrix Real) Any Any -> (Values (Matrix Real) (Listof Index)))
+                             ((Matrix Number)         -> (Values (Matrix Number) (Listof Index)))
+                             ((Matrix Number) Any     -> (Values (Matrix Number) (Listof Index)))
+                             ((Matrix Number) Any Any -> (Values (Matrix Number) (Listof Index)))))
+(define (matrix-gauss-elim M [jordan? #f] [unitize-pivot-row? #f])
   (define-values (m n) (matrix-shape M))
   (define rows (matrix->vector* M))
   (let loop ([#{i : Nonnegative-Fixnum} 0]
              [#{j : Nonnegative-Fixnum} 0]
-             [#{without-pivot : (Listof Index)}  '()])
+             [#{without-pivot : (Listof Index)} empty])
     (cond
-      [(and (i . fx< . m) (j . fx< . n))
-       ;; Find the row with the pivot value
-       (define p (cond [partial-pivoting?  (find-partial-pivot rows m i j)]
-                       [else  (find-pivot rows m i j)]))
-       (define pivot (if p (unsafe-vector2d-ref rows p j) 0))
-       (cond
-         [(or (not p) (zero? pivot))  ; didn't find pivot?
-          (loop i (fx+ j 1) (cons j without-pivot))]
-         [else 
-          (vector-swap! rows i p)  ; swap pivot row with current
-          (let ([pivot  (cond [reduced?  (vector-scale! (unsafe-vector-ref rows i) (/ pivot))
-                                         (/ pivot pivot)]
-                              [else  pivot])])
-            ;; Remove elements below pivot by scaling and adding the pivot's row to each row below
-            (let l-loop ([#{l : Nonnegative-Fixnum} (fx+ i 1)])
-              (cond [(l . fx< . m)
-                     (define x_lj (unsafe-vector2d-ref rows l j))
-                     (unless (zero? x_lj)
-                       (vector-scaled-add! (unsafe-vector-ref rows l)
-                                           (unsafe-vector-ref rows i)
-                                           (- (/ x_lj pivot))))
-                     (l-loop (fx+ l 1))]
-                    [else
-                     (loop (fx+ i 1) (fx+ j 1) without-pivot)])))])]
-      [else
+      [(j . fx>= . n)
        (values (vector*->matrix rows)
-               (reverse without-pivot))])))
+               (reverse without-pivot))]
+      [(i . fx>= . m)
+       (values (vector*->matrix rows)
+               ;; None of the rest of the columns can have pivots
+               (let loop ([#{j : Nonnegative-Fixnum} j] [without-pivot without-pivot])
+                 (cond [(j . fx< . n)  (loop (fx+ j 1) (cons j without-pivot))]
+                       [else  (reverse without-pivot)])))]
+      [else
+       (define-values (p pivot) (find-partial-pivot rows m i j))
+       (cond
+         [(zero? pivot)  (loop i (fx+ j 1) (cons j without-pivot))]
+         [else
+          ;; Swap pivot row with current
+          (vector-swap! rows i p)
+          ;; Possibly unitize the new current row
+          (let ([pivot  (if unitize-pivot-row?
+                            (begin (vector-scale! (unsafe-vector-ref rows i) (/ pivot))
+                                   1)
+                            pivot)])
+            (elim-rows! rows m i j pivot (if jordan? 0 (fx+ i 1)))
+            (loop (fx+ i 1) (fx+ j 1) without-pivot))])])))
 
-(: matrix-rank : (Matrix Number) -> Integer)
+;; ===================================================================================================
+;; Simple functions derived from Gaussian elimination
+
+(: matrix-row-echelon
+   (case-> ((Matrix Real)         -> (Matrix Real))
+           ((Matrix Real) Any     -> (Matrix Real))
+           ((Matrix Real) Any Any -> (Matrix Real))
+           ((Matrix Number)         -> (Matrix Number))
+           ((Matrix Number) Any     -> (Matrix Number))
+           ((Matrix Number) Any Any -> (Matrix Number))))
+(define (matrix-row-echelon M [jordan? #f] [unitize-pivot-row? jordan?])
+  (let-values ([(M _) (matrix-gauss-elim M jordan? unitize-pivot-row?)])
+    M))
+
+(: matrix-rank : (Matrix Number) -> Index)
 (define (matrix-rank M)
   ; TODO: Use QR or SVD instead for inexact matrices
   ; See answer: http://scicomp.stackexchange.com/questions/1861/understanding-how-numpy-does-svd
   ; rank = dimension of column space = dimension of row space
-  (define-values (m n) (matrix-shape M))
-  (define-values (_ cols-without-pivot) (matrix-gauss-eliminate M))
-  (- n (length cols-without-pivot)))
+  (define n (matrix-num-cols M))
+  (define-values (_ cols-without-pivot) (matrix-gauss-elim M))
+  (assert (- n (length cols-without-pivot)) index?))
 
-(: matrix-nullity : (Matrix Number) -> Integer)
+(: matrix-nullity : (Matrix Number) -> Index)
 (define (matrix-nullity M)
   ; nullity = dimension of null space
-  (define-values (m n) (matrix-shape M))
-  (define-values (_ cols-without-pivot) (matrix-gauss-eliminate M))
+  (define-values (_ cols-without-pivot)
+    (matrix-gauss-elim (ensure-matrix 'matrix-nullity M)))
   (length cols-without-pivot))
 
-(: matrix-determinant : (Matrix Number) -> Number)
+(: maybe-cons-submatrix (All (A) ((Matrix A) Nonnegative-Fixnum Nonnegative-Fixnum (Listof (Matrix A))
+                                             -> (Listof (Matrix A)))))
+(define (maybe-cons-submatrix M j0 j1 Bs)
+  (cond [(= j0 j1)  Bs]
+        [else  (cons (submatrix M (::) (:: j0 j1)) Bs)]))
+
+(: matrix-column-space (case-> ((Matrix Real) -> (Array Real))
+                               ((Matrix Number) -> (Array Number))))
+(define (matrix-column-space M)
+  (define n (matrix-num-cols M))
+  (define-values (_ wps) (matrix-gauss-elim M))
+  (cond [(empty? wps)  M]
+        [(= (length wps) n)  (make-array (vector 0 n) 0)]
+        [else
+         (define next-j (first wps))
+         (define Bs (maybe-cons-submatrix M 0 next-j empty))
+         (let loop ([#{j : Index} next-j] [wps (rest wps)] [Bs Bs])
+           (cond [(empty? wps)
+                  (matrix-augment (reverse (maybe-cons-submatrix M (fx+ j 1) n Bs)))]
+                 [else
+                  (define next-j (first wps))
+                  (loop next-j (rest wps) (maybe-cons-submatrix M (fx+ j 1) next-j Bs))]))]))
+
+;; ===================================================================================================
+;; Determinant
+
+(: matrix-determinant (case-> ((Matrix Real) -> Real)
+                              ((Matrix Number) -> Number)))
 (define (matrix-determinant M)
-  (define-values (m n) (matrix-shape M))
+  (define m (square-matrix-size M))
   (cond
-    [(= m 1) (matrix-ref M 0 0)]
-    [(= m 2) (let ([a (matrix-ref M 0 0)]
-                   [b (matrix-ref M 0 1)]
-                   [c (matrix-ref M 1 0)]
-                   [d (matrix-ref M 1 1)])
-               (- (* a d) (* b c)))]
-    [(= m 3) (let ([a (matrix-ref M 0 0)]
-                   [b (matrix-ref M 0 1)]
-                   [c (matrix-ref M 0 2)]
-                   [d (matrix-ref M 1 0)]
-                   [e (matrix-ref M 1 1)]
-                   [f (matrix-ref M 1 2)]
-                   [g (matrix-ref M 2 0)]
-                   [h (matrix-ref M 2 1)]
-                   [i (matrix-ref M 2 2)])
-               (+ (*    a  (- (* e i) (* f h)))
-                  (* (- b) (- (* d i) (* f g)))
-                  (*    c  (- (* d h) (* e g)))))]
-    [else           
-     (let-values ([(M _) (matrix-gauss-eliminate M #f #f)])
-       ; TODO: #f #f turns off partial pivoting
-       #; (for/product: : Number ([i (in-range 0 m)])
-            (matrix-ref M i i))
-       (let ()
-         (define: product : Number 1)
-         (for: ([i : Integer (in-range 0 m 1)])
-           (set! product (* product (matrix-ref M i i))))
-         product))]))
+    [(= m 1)  (matrix-ref M 0 0)]
+    [(= m 2)  (match-define (vector a b c d)
+                (mutable-array-data (array->mutable-array M)))
+              (- (* a d) (* b c))]
+    [(= m 3)  (match-define (vector a b c d e f g h i)
+                (mutable-array-data (array->mutable-array M)))
+              (+ (*    a  (- (* e i) (* f h)))
+                 (* (- b) (- (* d i) (* f g)))
+                 (*    c  (- (* d h) (* e g))))]
+    [else
+     (matrix-determinant/row-reduction M)]))
 
-(: matrix-column-space : (Matrix Number) -> (Listof (Matrix Number)))
-; Returns
-;  1) a list of column vectors spanning the column space
-;  2) a list of column vectors spanning the null space
-(define (matrix-column-space M)  
-  (define-values (m n) (matrix-shape M))
-  (: M1 (Matrix Number))
-  (: cols-without-pivot (Listof Integer))
-  (define-values (M1 cols-without-pivot) (matrix-gauss-eliminate M #t))
-  (set! M1 (array->mutable-array M1))  
-  (define: column-space : (Listof (Matrix Number))
-    (for/list:
-        ([i : Index n]
-         #:when (not (member i cols-without-pivot)))
-      (matrix-col M1 i)))
-  column-space)
-
-(: matrix-row-echelon-form : 
-   (case-> ((Matrix Number) Boolean -> (Matrix Number))
-           ((Matrix Number) Boolean -> (Matrix Number))
-           ((Matrix Number)         -> (Matrix Number))))
-(define (matrix-row-echelon-form M [unitize-pivot-row? #f])
-  (let-values ([(M wp) (matrix-gauss-eliminate M unitize-pivot-row?)])
-    M))
-
-(: matrix-gauss-jordan-eliminate : 
-   (case-> ((Matrix Number) Boolean Boolean -> (Values (Matrix Number) (Listof Integer)))
-           ((Matrix Number) Boolean         -> (Values (Matrix Number) (Listof Integer)))
-           ((Matrix Number)                 -> (Values (Matrix Number) (Listof Integer)))))
-(define (matrix-gauss-jordan-eliminate M [unitize-pivot-row? #f] [partial-pivoting? #t])
-  (define-values (m n) (matrix-shape M))
-  (: loop : (Integer Integer (Matrix Number) Integer (Listof Integer) 
-                     -> (Values (Matrix Number) (Listof Integer))))
-  (define (loop i j ; i from 0 to m
-                M
-                k   ; count rows without pivot
-                without-pivot) 
+(: matrix-determinant/row-reduction (case-> ((Matrix Real) -> Real)
+                                            ((Matrix Number) -> Number)))
+(define (matrix-determinant/row-reduction M)
+  (define m (square-matrix-size M))
+  (define rows (matrix->vector* M))
+  (let loop ([#{i : Nonnegative-Fixnum} 0] [#{sign : Real} 1])
     (cond
-      [(or (= i m) (= j n)) (values M without-pivot)]
-      [else
-       ; find row to become pivot
-       (define p
-         (if partial-pivoting?
-             ; find element with maximal absolute value
-             (let: max-loop : (U False Integer) 
-               ([l : Integer i] ; i<=l<m
-                [max-current : Real -inf.0]
-                [max-index : Integer i])
-               (cond
-                 [(= l m) max-index]
-                 [else 
-                  (let ([v (magnitude (matrix-ref M l j))])
-                    (if (> (magnitude (matrix-ref M l j)) max-current)
-                        (max-loop (+ l 1) v l)
-                        (max-loop (+ l 1) max-current max-index)))]))
-             ; find non-zero element in column
-             (let: first-loop : (U False Integer) 
-               ([l : Integer i]) ; i<=l<m
-               (cond
-                 [(= l m) #f]
-                 [(not (zero? (matrix-ref M l j))) l]
-                 [else (first-loop (+ l 1))]))))
+      [(i . fx< . m)
+       (define-values (p pivot) (find-partial-pivot rows m i i))
        (cond
-         [(eq? p #f)
-          ; no pivot found - this implies the matrix is singular (not invertible)
-          (loop i (+ j 1) M (+ k 1) (cons j without-pivot))]
-         [else 
-          ; swap if neccessary
-          (let* ([M (if (= i p) M (matrix-swap-rows M i p))]
-                 ; now we now (i,j) is a pivot
-                 [M ; maybe scale row
-                  (if unitize-pivot-row?
-                      (let ([pivot (matrix-ref M i j)])
-                        (if (zero? pivot)
-                            M
-                            (matrix-scale-row M i (/ pivot))))
-                      M)])
-            (let ([pivot (matrix-ref M i j)])
-              ; remove elements above and below pivot
-              (let l-loop ([l 0] [M M])
-                (cond
-                  [(= l m) (loop (+ i 1) (+ j 1) M k without-pivot)]
-                  [(= l i) (l-loop (+ l 1) M)]
-                  [else
-                   (let ([x_lj (matrix-ref M l j)])
-                     (l-loop (+ l 1)
-                             (if (zero? x_lj)
-                                 M
-                                 (matrix-add-scaled-row M l (- (/ x_lj pivot)) i))))]))))])]))
-  (let-values ([(M without) (loop 0 0 M 0 '())])
-    (values M without)))
+         [(zero? pivot)  0]  ; no pivot means non-invertible matrix
+         [else
+          (vector-swap! rows i p)                  ; negates determinant if i != p
+          (elim-rows! rows m i i pivot (fx+ i 1))  ; doesn't change the determinant
+          (loop (fx+ i 1) (if (= i p) sign (* -1 sign)))])]
+      [else
+       (define prod (unsafe-vector2d-ref rows 0 0))
+       (let loop ([#{i : Nonnegative-Fixnum} 1] [prod prod])
+         (cond [(i . fx< . m)
+                (loop (fx+ i 1) (* prod (unsafe-vector2d-ref rows i i)))]
+               [else  (* prod sign)]))])))
 
-(: matrix-reduced-row-echelon-form : 
-   (case-> ((Matrix Number) Boolean -> (Matrix Number))
-           ((Matrix Number) Boolean -> (Matrix Number))
-           ((Matrix Number)         -> (Matrix Number))))
-(define (matrix-reduced-row-echelon-form M [unitize-pivot-row? #f])
-  (let-values ([(M wp) (matrix-gauss-jordan-eliminate M unitize-pivot-row?)])
-    M))
+;; ===================================================================================================
+;; Inversion and solving linear systems
 
-(: matrix-inverse : (Matrix Number) -> (Matrix Number))
-(define (matrix-inverse M)
-  (define-values (m n) (matrix-shape M))
-  (unless (= m n) (error 'matrix-inverse "matrix not square"))
-  (let ([MI (matrix-augment (list M (identity-matrix m)))])
-    (define 2m (* 2 m))
-    (if (index? 2m)
-        (submatrix (matrix-reduced-row-echelon-form MI #t) 
-                   (in-range 0 m) (in-range m 2m))
-        (error 'matrix-inverse "internal error"))))
+(: matrix-invertible? ((Matrix Number) -> Boolean))
+(define (matrix-invertible? M)
+  (not (zero? (matrix-determinant M))))
 
-(: matrix-solve : (Matrix Number) (Matrix Number) -> (Matrix Number))
-;    Return a column-vector x such that Mx = b.
-;    If no such vector exists return #f.
-(define (matrix-solve M b)
-  (define-values (m n) (matrix-shape M))
-  (define-values (s t) (matrix-shape b))
-  (define m+1 (+ m 1))
-  (cond
-    [(not (= t 1)) (error 'matrix-solve "expected column vector (i.e. r x 1 - matrix), got: ~a " b)]
-    [(not (= m s)) (error 'matrix-solve "expected column vector with same number of rows as the matrix")]
-    [(index? m+1)
-     (submatrix
-      (matrix-reduced-row-echelon-form 
-       (matrix-augment (list M b)) #t)
-      (in-range 0 m) (in-range m m+1))]
-    [else (error 'matrix-solve "internatl error")]))
+(: make-invertible-fail (Symbol (Matrix Any) -> (-> Nothing)))
+(define ((make-invertible-fail name M))
+  (raise-argument-error name "matrix-invertible?" M))
 
-(: matrix-solve-many : (Matrix Number) (Listof (Matrix Number)) -> (Matrix Number))
-(define (matrix-solve-many M bs)
-  (define-values (m n) (matrix-shape M))
-  (define-values (s t) (matrix-shape (car bs)))
-  (define k (length bs))
-  (define m+1 (+ m 1))
-  (define m+k (+ m k))
-  (cond
-    [(not (= t 1)) (error 'matrix-solve-many "expected column vector (i.e. r x 1 - matrix), got: ~a " (car bs))]
-    [(not (= m s)) (error 'matrix-solve-many "expected column vectors with same number of rows as the matrix")]
-    [(and (index? m+1) (index? m+k)) 
-     (define bs-as-matrix (matrix-augment bs))     
-     (define MB (matrix-augment (list M bs-as-matrix)))
-     (define reduced-MB (matrix-reduced-row-echelon-form MB #t))
-     (submatrix reduced-MB 
-                (in-range 0 m+k)
-                (in-range m m+1))]
-    [else (error 'matrix-solve-many "internal error")]))
+(: matrix-inverse (All (A) (case-> ((Matrix Real)        -> (Matrix Real))
+                                   ((Matrix Real) (-> A) -> (U A (Matrix Real)))
+                                   ((Matrix Number)        -> (Matrix Number))
+                                   ((Matrix Number) (-> A) -> (U A (Matrix Number))))))
+(define matrix-inverse
+  (case-lambda
+    [(M)  (matrix-inverse M (make-invertible-fail 'matrix-inverse M))]
+    [(M fail)
+     (define m (square-matrix-size M))
+     (define I (identity-matrix m))
+     (define-values (IM^-1 wps) (matrix-gauss-elim (matrix-augment (list M I)) #t #t))
+     (cond [(and (not (empty? wps)) (= (first wps) m))
+            (submatrix IM^-1 (::) (:: m #f))]
+           [else  (fail)])]))
 
+(: matrix-solve (All (A) (case->
+                          ((Matrix Real) (Matrix Real)        -> (Matrix Real))
+                          ((Matrix Real) (Matrix Real) (-> A) -> (U A (Matrix Real)))
+                          ((Matrix Number) (Matrix Number)        -> (Matrix Number))
+                          ((Matrix Number) (Matrix Number) (-> A) -> (U A (Matrix Number))))))
+(define matrix-solve
+  (case-lambda
+    [(M B)  (matrix-solve M B (make-invertible-fail 'matrix-solve M))]
+    [(M B fail)
+     (define m (square-matrix-size M))
+     (define-values (s t) (matrix-shape B))
+     (cond [(= m s)
+            (define-values (IX wps) (matrix-gauss-elim (matrix-augment (list M B)) #t #t))
+            (cond [(and (not (empty? wps)) (= (first wps) m))
+                   (submatrix IX (::) (:: m #f))]
+                  [else  (fail)])]
+           [else
+            (error 'matrix-solve
+                   "matrices must have the same number of rows; given ~e and ~e"
+                   M B)])]))
 
-;;; LU Factorization
-; Not all matrices can be LU-factored.
-; If Gauss-elimination can be done without any row swaps,
-; a LU-factorization is possible.
+;; ===================================================================================================
+;; LU Factorization
+
+;; An LU factorization exists iff Gaussian elimination can be done without row swaps.
 
 (: matrix-lu : 
    (Matrix Number) -> (U False (List (Matrix Number) (Matrix Number))))
@@ -625,6 +574,9 @@
       (take (sort (loop vs '() 0) norm>) r)
       (error 'extend-span-to-basis "expected index as second argument, got ~a" r)))
 
+;; ===================================================================================================
+;; QR decomposition
+
 (: matrix-qr : (Matrix Number) -> (Values (Matrix Number) (Matrix Number)))
 (define (matrix-qr M)
   ; compute the QR-facorization
diff --git a/collects/math/private/matrix/matrix-types.rkt b/collects/math/private/matrix/matrix-types.rkt
index 74925601fa..0eaf70bf29 100644
--- a/collects/math/private/matrix/matrix-types.rkt
+++ b/collects/math/private/matrix/matrix-types.rkt
@@ -71,9 +71,9 @@
 (: matrix-num-rows (All (A) ((Array A) -> Index)))
 (define (matrix-num-rows a)
   (cond [(matrix? a)  (vector-ref (array-shape a) 0)]
-        [else  (raise-argument-error 'matrix-col-length "matrix?" a)]))
+        [else  (raise-argument-error 'matrix-num-rows "matrix?" a)]))
 
 (: matrix-num-cols (All (A) ((Array A) -> Index)))
 (define (matrix-num-cols a)
   (cond [(matrix? a)  (vector-ref (array-shape a) 1)]
-        [else  (raise-argument-error 'matrix-row-length "matrix?" a)]))
+        [else  (raise-argument-error 'matrix-num-cols "matrix?" a)]))
diff --git a/collects/math/tests/matrix-tests.rkt b/collects/math/tests/matrix-tests.rkt
index f8d7716574..24febab056 100644
--- a/collects/math/tests/matrix-tests.rkt
+++ b/collects/math/tests/matrix-tests.rkt
@@ -20,7 +20,7 @@
         (make-array #(0 1) 0)
         (make-array #(0 0) 0)
         (make-array #(1 1 1) 0)))
-#|
+
 ;; ===================================================================================================
 ;; Literal syntax
 
@@ -666,32 +666,136 @@
 (check-exn exn:fail:contract? (λ () (matrix-trace (col-matrix [1 2 3]))))
 (for: ([a  (in-list nonmatrices)])
   (check-exn exn:fail:contract? (λ () (matrix-trace a))))
-|#
+
 ;; ===================================================================================================
 ;; Gaussian elimination
 
-(: gauss-eliminate : (Matrix Number) Boolean Boolean -> (Matrix Number))
-(define (gauss-eliminate M reduce? partial-pivot?)
-  (let-values ([(M wp) (matrix-gauss-eliminate M reduce? partial-pivot?)])
-    M))
+(check-equal? (matrix-row-echelon (matrix [[2 4] [3 4]]) #f #f)
+              (matrix [[3 4] [0 4/3]]))
 
-(check-equal? (gauss-eliminate (matrix [[1 2] [3 4]]) #f #f)
-              (matrix [[1 2] [0 -2]]))
+(check-equal? (matrix-row-echelon (matrix [[2 4] [3 4]]) #f #t)
+              (matrix [[1 4/3] [0 1]]))
 
-(check-equal? (gauss-eliminate (matrix [[2 4] [3 4]]) #t #f)
-              (matrix [[1 2] [0 1]]))
-
-(check-equal? (gauss-eliminate (matrix [[2. 4.] [3. 4.]]) #t #t)
-              (matrix [[1. 1.3333333333333333] [0. 1.]]))
-
-(check-equal? (gauss-eliminate (matrix [[1 4] [2 4]]) #t #t)
-              (matrix [[1 2] [0 1]]))
-
-(check-equal? (gauss-eliminate (matrix [[1 2] [2 4]]) #f #t)
+(check-equal? (matrix-row-echelon (matrix [[1 2] [2 4]]) #f #f)
               (matrix [[2 4] [0 0]]))
 
+(check-equal? (matrix-row-echelon (matrix [[1 4] [2 4]]) #f #t)
+              (matrix [[1 2] [0 1]]))
+
+(check-equal? (matrix-row-echelon (matrix [[ 2  1 -1   8]
+                                           [-3 -1  2 -11]
+                                           [-2  1  2  -3]])
+                                  #f #t)
+              (matrix [[1 1/3 -2/3 11/3]
+                       [0  1   2/5 13/5]
+                       [0  0    1   -1]]))
+
+(check-equal? (matrix-row-echelon (matrix [[ 2  1 -1   8]
+                                           [-3 -1  2 -11]
+                                           [-2  1  2  -3]])
+                                  #t)
+              (matrix [[1 0 0  2]
+                       [0 1 0  3]
+                       [0 0 1 -1]]))
+
 (for: ([a  (in-list nonmatrices)])
-  (check-exn exn:fail:contract? (λ () (gauss-eliminate a #f #f))))
+  (check-exn exn:fail:contract? (λ () (matrix-row-echelon a))))
+
+(check-equal? (matrix-rank (matrix [[0 0] [0 0]])) 0)
+(check-equal? (matrix-rank (matrix [[1 0] [0 0]])) 1)
+(check-equal? (matrix-rank (matrix [[1 0] [0 3]])) 2)
+(check-equal? (matrix-rank (matrix [[1 2] [2 4]])) 1)
+(check-equal? (matrix-rank (matrix [[1 2] [3 4]])) 2)
+(check-equal? (matrix-rank (matrix [[1 2 3]])) 1)
+(check-equal? (matrix-rank (matrix [[1 2 3] [2 3 5]])) 2)
+(check-equal? (matrix-rank (matrix [[1 2 3] [2 3 5] [3 4 7]])) 2)
+(check-equal? (matrix-rank (matrix [[1 2 3] [2 3 5] [3 4 7] [4 5 9]])) 2)
+(check-equal? (matrix-rank (matrix [[1 2 3 5] [2 3 5 8]])) 2)
+(check-equal? (matrix-rank (matrix [[1 5 2 3] [2 8 3 5]])) 2)
+(for: ([a  (in-list nonmatrices)])
+  (check-exn exn:fail:contract? (λ () (matrix-rank a))))
+
+(check-equal? (matrix-nullity (matrix [[0 0] [0 0]])) 2)
+(check-equal? (matrix-nullity (matrix [[1 0] [0 0]])) 1)
+(check-equal? (matrix-nullity (matrix [[1 0] [0 3]])) 0)
+(check-equal? (matrix-nullity (matrix [[1 2] [2 4]])) 1)
+(check-equal? (matrix-nullity (matrix [[1 2] [3 4]])) 0)
+(check-equal? (matrix-nullity (matrix [[1 2 3]])) 2)
+(check-equal? (matrix-nullity (matrix [[1 2 3] [2 3 5]])) 1)
+(check-equal? (matrix-nullity (matrix [[1 2 3] [2 3 5] [3 4 7]])) 1)
+(check-equal? (matrix-nullity (matrix [[1 2 3] [2 3 5] [3 4 7] [4 5 9]])) 1)
+(check-equal? (matrix-nullity (matrix [[1 2 3 5] [2 3 5 8]])) 2)
+(check-equal? (matrix-nullity (matrix [[1 5 2 3] [2 8 3 5]])) 2)
+(for: ([a  (in-list nonmatrices)])
+  (check-exn exn:fail:contract? (λ () (matrix-nullity a))))
+
+;; ===================================================================================================
+;; Determinant
+
+(check-equal? (matrix-determinant (matrix [[3]])) 3)
+(check-equal? (matrix-determinant (matrix [[1 2] [3 4]])) (- (* 1 4) (* 2 3)))
+(check-equal? (matrix-determinant (matrix [[1 2 3] [4  5 6] [7 8 9]])) 0)
+(check-equal? (matrix-determinant (matrix [[1 2 3] [4 -5 6] [7 8 9]])) 120)
+(check-equal? (matrix-determinant (matrix [[1 2 3 4]
+                                           [-5 6 7 8]
+                                           [9 10 -11 12]
+                                           [13 14 15 16]]))
+              5280)
+
+(for: ([_  (in-range 100)])
+  (define a (array- (random-matrix 3 3 7) (array 3)))
+  (check-equal? (matrix-determinant/row-reduction a)
+                (matrix-determinant a)))
+
+(check-exn exn:fail:contract? (λ () (matrix-determinant (matrix [[1 2 3] [4 5 6]]))))
+(check-exn exn:fail:contract? (λ () (matrix-determinant (matrix [[1 4] [2 5] [3 6]]))))
+(for: ([a  (in-list nonmatrices)])
+  (check-exn exn:fail:contract? (λ () (matrix-determinant a))))
+
+;; ===================================================================================================
+;; Solving linear systems
+
+(for: ([_  (in-range 100)])
+  (define M (array- (random-matrix 3 3 7) (array 3)))
+  (define B (array- (random-matrix 3 (+ 1 (random 10)) 7) (array 3)))
+  (cond [(matrix-invertible? M)
+         (define X (matrix-solve M B))
+         (check-equal? (matrix* M X) B (format "M = ~a  B = ~a" M B))]
+        [else
+         (check-false (matrix-solve M B (λ () #f))
+                      (format "M = ~a  B = ~a" M B))]))
+
+(check-exn exn:fail? (λ () (matrix-solve (random-matrix 3 4) (random-matrix 3 1))))
+(check-exn exn:fail? (λ () (matrix-solve (random-matrix 4 3) (random-matrix 4 1))))
+
+(check-exn exn:fail:contract? (λ () (matrix-solve (random-matrix 3 4) (random-matrix 4 1))))
+(check-exn exn:fail:contract? (λ () (matrix-solve (random-matrix 4 3) (random-matrix 3 1))))
+
+(for: ([a  (in-list nonmatrices)])
+  (check-exn exn:fail:contract? (λ () (matrix-solve a (matrix [[1]]))))
+  (check-exn exn:fail:contract? (λ () (matrix-solve (matrix [[1]]) a))))
+
+;; ===================================================================================================
+;; Inversion
+
+(for: ([_  (in-range 100)])
+  (define a (array- (random-matrix 3 3 7) (array 3)))
+  (cond [(matrix-invertible? a)
+         (check-equal? (matrix* a (matrix-inverse a))
+                       (identity-matrix 3)
+                       (format "~a" a))
+         (check-equal? (matrix* (matrix-inverse a) a)
+                       (identity-matrix 3)
+                       (format "~a" a))]
+        [else
+         (check-false (matrix-inverse a (λ () #f))
+                      (format "~a" a))]))
+
+(check-exn exn:fail:contract? (λ () (matrix-inverse (random-matrix 3 4))))
+(check-exn exn:fail:contract? (λ () (matrix-inverse (random-matrix 4 3))))
+
+(for: ([a  (in-list nonmatrices)])
+  (check-exn exn:fail:contract? (λ () (matrix-inverse a))))
 
 #|
 ;; ===================================================================================================
@@ -773,37 +877,6 @@
                       4 4 ((inst vector Number)
                            2 4 9 11 0 0.0 2.23606797749979 2.23606797749979 
                            0 0 0.0 4.440892098500626e-16 0 0 0 0.0))))))
-    (list 'matrix-solve
-          (let* ([M (list*->matrix '[[1 5] [2 3]])] 
-                 [b (list*->matrix '[[5] [5]])])
-            (equal? (matrix* M (matrix-solve M b)) b)))
-    (list 'matrix-inverse
-          (equal? (let ([M (list*->matrix '[[1 2] [3 4]])]) (matrix* M (matrix-inverse M)))
-                  (identity-matrix 2))
-          (equal? (let ([M (list*->matrix '[[1 2] [3 4]])]) (matrix* (matrix-inverse M) M))
-                  (identity-matrix 2)))
-    (list 'matrix-determinant
-          (equal? (matrix-determinant (list*->matrix '[[3]])) 3)
-          (equal? (matrix-determinant (list*->matrix '[[1 2] [3 4]])) (- (* 1 4) (* 2 3)))
-          (equal? (matrix-determinant (list*->matrix '[[1 2 3] [4  5 6] [7 8 9]])) 0)
-          (equal? (matrix-determinant (list*->matrix '[[1 2 3] [4 -5 6] [7 8 9]])) 120)
-          (equal? (matrix-determinant (list*->matrix '[[1 2 3 4]
-                                                       [-5 6 7 8]
-                                                       [9 10 -11 12]
-                                                       [13 14 15 16]]))
-                  5280))
-    (list 
-     'matrix-scale-row
-     (equal? (matrix-scale-row (identity-matrix 3) 0 2)
-             (list*->array '[[2 0 0] [0 1 0] [0 0 1]] real? )))
-    (list
-     'matrix-swap-rows
-     (equal? (matrix-swap-rows (list*->array '[[1 2 3] [4 5 6] [7 8 9]] real? ) 0 1)
-             (list*->array '[[4 5 6] [1 2 3] [7 8 9]] real? )))
-    (list
-     'matrix-add-scaled-row
-     (equal? (matrix-add-scaled-row (list*->array '[[1 2 3] [4 5 6] [7 8 9]] real? ) 0 2 1)
-             (list*->array '[[9 12 15] [4 5 6] [7 8 9]] real? )))
     (let ()
       (define M (list*->matrix '[[1  1  0  3]
                                  [2  1 -1  1]
@@ -828,39 +901,6 @@
                           [0  0  3  13]
                           [0  0  0 -13]]))
              (equal? (matrix* L V) M)))))
-    (list 
-     'matrix-rank
-     (equal? (matrix-rank (list*->matrix '[[0 0] [0 0]])) 0)
-     (equal? (matrix-rank (list*->matrix '[[1 0] [0 0]])) 1)
-     (equal? (matrix-rank (list*->matrix '[[1 0] [0 3]])) 2)
-     (equal? (matrix-rank (list*->matrix '[[1 2] [2 4]])) 1)
-     (equal? (matrix-rank (list*->matrix '[[1 2] [3 4]])) 2))
-    (list 
-     'matrix-nullity
-     (equal? (matrix-nullity (list*->matrix '[[0 0] [0 0]])) 2)
-     (equal? (matrix-nullity (list*->matrix '[[1 0] [0 0]])) 1)
-     (equal? (matrix-nullity (list*->matrix '[[1 0] [0 3]])) 0)
-     (equal? (matrix-nullity (list*->matrix '[[1 2] [2 4]])) 1)
-     (equal? (matrix-nullity (list*->matrix '[[1 2] [3 4]])) 0))
-    #;
-    (let ()
-      (define-values (c1 n1) 
-        (matrix-column+null-space (list*rix '[[0 0] [0 0]])))
-      (define-values (c2 n2) 
-        (matrix-column+null-space (list*->matrix '[[1 2] [2 4]])))
-      (define-values (c3 n3) 
-        (matrix-column+null-space (list*atrix '[[1 2] [2 5]])))
-      (list 
-       'matrix-column+null-space
-       (equal? c1 '())
-       (equal? n1 (list (list*->matrix '[[0] [0]])
-                        (list*->matrix '[[0] [0]])))
-       (equal? c2 (list (list*->matrix '[[1] [2]])))
-       ;(equal? n2 '([0 0]))
-       (equal? c3 (list (list*->matrix '[[1] [2]])
-                        (list*->matrix '[[2] [5]])))
-       (equal? n3 '()))))
-  
   #;
   (begin
     "matrix-2d.rkt"