diff --git a/collects/math/private/matrix/matrix-basic.rkt b/collects/math/private/matrix/matrix-basic.rkt
index 5ed6a2f49f..2556544c96 100644
--- a/collects/math/private/matrix/matrix-basic.rkt
+++ b/collects/math/private/matrix/matrix-basic.rkt
@@ -70,7 +70,8 @@
         [else
          (unsafe-array-ref a ((inst vector Index) i j))]))
 
-(: submatrix (All (A) (Matrix A) Slice-Spec Slice-Spec -> (Matrix A)))
+(: submatrix (All (A) (Matrix A) (U Slice (Sequenceof Integer)) (U Slice (Sequenceof Integer))
+                  -> (Array A)))
 (define (submatrix a row-range col-range)
   (array-slice-ref (ensure-matrix 'submatrix a) (list row-range col-range)))
 
@@ -120,11 +121,11 @@
 
 (: matrix-diagonal (All (A) ((Matrix A) -> (Array A))))
 (define (matrix-diagonal a)
-  (define m (square-matrix-size a))
+  (define-values (m n) (matrix-shape a))
   (define proc (unsafe-array-proc a))
   (array-default-strict
    (unsafe-build-array
-    ((inst vector Index) m)
+    ((inst vector Index) (fxmin m n))
     (λ: ([js : Indexes])
       (define i (unsafe-vector-ref js 0))
       (proc ((inst vector Index) i i))))))
@@ -305,8 +306,11 @@
 (: matrix-trace (case-> ((Matrix Real) -> Real)
                         ((Matrix Number) -> Number)))
 (define (matrix-trace a)
-  (parameterize ([array-strictness #f])
-    (array-all-sum (matrix-diagonal a))))
+  (cond [(square-matrix? a)
+         (parameterize ([array-strictness #f])
+           (array-all-sum (matrix-diagonal a)))]
+        [else
+         (raise-argument-error 'matrix-trace "square-matrix?" a)]))
 
 ;; ===================================================================================================
 ;; Row/column operations
diff --git a/collects/math/private/matrix/matrix-solve.rkt b/collects/math/private/matrix/matrix-solve.rkt
index f5a7d66542..e92fdaa120 100644
--- a/collects/math/private/matrix/matrix-solve.rkt
+++ b/collects/math/private/matrix/matrix-solve.rkt
@@ -69,7 +69,8 @@
 
 (: matrix-invertible? ((Matrix Number) -> Boolean))
 (define (matrix-invertible? M)
-  (not (zero? (matrix-determinant M))))
+  (and (square-matrix? M)
+       (not (zero? (matrix-determinant M)))))
 
 (: matrix-inverse (All (A) (case-> ((Matrix Real)        -> (Matrix Real))
                                    ((Matrix Real) (-> A) -> (U A (Matrix Real)))
diff --git a/collects/math/scribblings/math-matrix.scrbl b/collects/math/scribblings/math-matrix.scrbl
index 94039bcee8..f891d4b8a1 100644
--- a/collects/math/scribblings/math-matrix.scrbl
+++ b/collects/math/scribblings/math-matrix.scrbl
@@ -6,7 +6,7 @@
                      racket/list
                      math plot
                      (only-in typed/racket/base
-                              ann inst : λ: define: make-predicate ->
+                              ann inst : λ: define: make-predicate -> case-> Nothing
                               Flonum Real Boolean Any Integer Index Natural Exact-Positive-Integer
                               Nonnegative-Real Sequenceof Fixnum Values Number Float-Complex
                               All U List Vector Listof Vectorof Struct FlVector
@@ -41,9 +41,16 @@ For now, if you need speed, use the @racketmodname[typed/racket] language.
 
 @defmodule[math/matrix]
 
-Documentation for this module is currently under construction.
-
-Intro topics: definitions, case-> types, non-strictness
+Like all of @racketmodname[math], @racketmodname[math/matrix] is a work in progress.
+Most of the basic algorithms are implemented, but some are still in planning.
+Possibly the most useful unimplemented algorithms are
+@itemlist[@item{LUP decomposition (currently, LU decomposition is implemented, in @racket[matrix-lu])}
+          @item{@racket[matrix-solve] for upper-triangular matrices}
+          @item{Singular value decomposition (SVD)}
+          @item{Eigendecomposition}
+          @item{Decomposition-based solvers}
+          @item{Pseudoinverse, least-squares fitting}]
+All of these are planned for the next Racket release.
 
 @local-table-of-contents[]
 
@@ -51,6 +58,76 @@ Intro topics: definitions, case-> types, non-strictness
 @;{==================================================================================================}
 
 
+@section[#:tag "matrix:intro"]{Introduction}
+
+From the point of view of the functions in @racketmodname[math/matrix], a @deftech{matrix} is an
+@racket[Array] with two axes and at least one entry, or an array for which @racket[matrix?] returns
+@racket[#t].
+
+Technically, a matrix's entries may be any type, and some fully polymorphic matrix functions such as
+@racket[matrix-row] and @racket[matrix-map] operate on any kind of matrix.
+Other functions, such as @racket[matrix+], require their input matrices to contain either
+@racket[Real] or @racket[Number] values.
+
+@subsection[#:tag "matrix:function-types"]{Function Types}
+
+The documentation for @racketmodname[math/matrix] functions use the type @racket[Matrix], a synonym
+of @racket[Array], when the function either requires that an argument is a @tech{matrix} or ensures
+that a return value is a matrix.
+
+Most functions that implement matrix algorithms are documented as accepting @racket[(Matrix Number)]
+values. This includes @racket[(Matrix Real)], which is a subtype. Most of these functions have a more
+precise type than is documented. For example, @racket[matrix-conjugate] actually has the type
+@racketblock[(case-> ((Matrix Real)   -> (Matrix Real))
+                     ((Matrix Number) -> (Matrix Number)))]
+even though it is documented as having the less precise type
+@racket[((Matrix Number) -> (Matrix Number))].
+
+Precise function types allow Typed Racket to prove more facts about @racketmodname[math/matrix]
+client programs. In particular, it is usually easy for it to prove that matrix expressions on real
+matrices return real matrices:
+@interaction[#:eval typed-eval
+                    (matrix-conjugate (matrix [[1 2 3] [4 5 6]]))]
+
+@subsection[#:tag "matrix:failure"]{Failure Arguments}
+
+In many matrix operations, such as inversion, failure is easy to detect during computation, but
+is just as expensive to detect ahead of time as the operation itself.
+In these cases, the functions implementing the operations accept an optional @deftech{failure thunk},
+or a zero-argument function that returns the result of the operation in case of failure.
+
+For example, the (simplified) type of @racket[matrix-inverse] is
+@racketblock[(All (F) (case-> ((Matrix Number)        -> (Matrix Number))
+                              ((Matrix Number) (-> F) -> (U F (Matrix Number)))))]
+Thus, if a failure thunk is given, the call site is required to check for return values of type
+@racket[F] explicitly.
+
+The default failure thunk, which raises an error, has type @racket[(-> Nothing)].
+For such failure thunks, @racket[(U F (Matrix Number))] is equivalent to @racket[(Matrix Number)],
+because @racket[Nothing] is part of every type. (In Racket, any expression may raise an error.)
+Thus, in this case, no explicit test for values of type @racket[F] is necessary (though of course they
+may be caught using @racket[with-handlers] or similar).
+
+@subsection[#:tag "matrix:broadcasting"]{Broadcasting}
+
+Pointwise matrix operations do not @tech{broadcast} their arguments when given matrices with
+different sizes:
+@interaction[#:eval typed-eval
+                    (matrix+ (identity-matrix 2) (matrix [[10]]))]
+If you need broadcasting, use array operations:
+@interaction[#:eval typed-eval
+                    (array+ (identity-matrix 2) (matrix [[10]]))]
+
+@subsection{Strictness}
+
+Functions exported by @racketmodname[math/matrix] return @tech{strict} or @tech{nonstrict} arrays
+based on the value of the @racket[array-strictness] parameter.
+See @secref{array:nonstrict} for details.
+
+
+@;{==================================================================================================}
+
+
 @section[#:tag "matrix:types"]{Types, Predicates and Accessors}
 
 @defform[(Matrix A)]{
@@ -59,7 +136,7 @@ Equivalent to @racket[(Array A)], but used for values @racket[M] for which @rack
 }
 
 @defproc[(matrix? [arr (Array A)]) Boolean]{
-Returns @racket[#t] when @racket[arr] is a @deftech{matrix}: a nonempty array with exactly two axes.
+Returns @racket[#t] when @racket[arr] is a @tech{matrix}: a nonempty array with exactly two axes.
 @examples[#:eval typed-eval
                  (matrix? (array 10))
                  (matrix? (array #[1 2 3]))
@@ -161,7 +238,9 @@ Returns a matrix with @racket[xs] along the diagonal and @racket[0] everywhere e
 The length of @racket[xs] must be positive.
 }
 
-@margin-note{@hyperlink["http://en.wikipedia.org/wiki/Block_matrix#Block_diagonal_matrices"]{Wikipedia: Block-diagonal matrices}}
+@define[block-diagonal-url]{http://en.wikipedia.org/wiki/Block_matrix#Block_diagonal_matrices}
+
+@margin-note*{@hyperlink[block-diagonal-url]{Wikipedia: Block-diagonal matrices}}
 @defproc[(block-diagonal-matrix [Xs (Listof (Matrix A))]) (Matrix (U A 0))]{
 Returns a matrix with matrices @racket[Xs] along the diagonal and @racket[0] everywhere else.
 The length of @racket[Xs] must be positive.
@@ -172,7 +251,9 @@ The length of @racket[Xs] must be positive.
                                               (row-matrix [4 5 6])))]
 }
 
-@margin-note{@hyperlink["http://en.wikipedia.org/wiki/Vandermonde_matrix"]{Wikipedia: Vandermonde matrix}}
+@define[vandermonde-url]{http://en.wikipedia.org/wiki/Vandermonde_matrix}
+
+@margin-note*{@hyperlink[vandermonde-url]{Wikipedia: Vandermonde matrix}}
 @defproc[(vandermonde-matrix [xs (Listof Number)] [n Integer]) (Matrix Number)]{
 Returns an @racket[m]×@racket[n] Vandermonde matrix, where @racket[m = (length xs)].
 @examples[#:eval typed-eval
@@ -313,11 +394,11 @@ where each entry in @racket[M] is multiplied with @racket[z].
                  (matrix-scale (matrix [[1 2] [3 4]]) 2)]
 }
 
-@defproc*[([(matrix-map [f (A -> R)] [arr0 (Matrix A)]) (Matrix R)]
-           [(matrix-map [f (A B Ts ... -> R)] [arr0 (Matrix A)] [arr1 (Matrix B)] [arrs (Matrix Ts)]
+@defproc*[([(matrix-map [f (A -> R)] [M (Matrix A)]) (Matrix R)]
+           [(matrix-map [f (A B Ts ... -> R)] [M0 (Matrix A)] [M1 (Matrix B)] [N (Matrix Ts)]
                         ...)
             (Matrix R)])]{
-Like @racket[array-map], but requires at least one array argument and never @tech{broadcasts}.
+Like @racket[array-map], but requires at least one matrix argument and never @tech{broadcasts}.
 @examples[#:eval untyped-eval
                  (matrix-map sqr (matrix [[1 2] [3 4]]))
                  (matrix-map  +  (matrix [[1 2] [3 4]])
@@ -344,47 +425,52 @@ tolerant to floating-point error.
 @defproc[(matrix-ref [M (Matrix A)] [i Integer] [j Integer]) A]{
 Returns the entry on row @racket[i] and column @racket[j].
 @examples[#:eval untyped-eval
-                 (define A (matrix ([1 2 3] [4 5 6])))
-                 (matrix-ref A 0 2)
-                 (matrix-ref A 1 2)]
-}
-
-@defthing[submatrix Procedure]{
-@;{ TODO
-(: submatrix (All (A) (Matrix A) Slice-Spec Slice-Spec -> (Matrix A)))
-(define (submatrix a row-range col-range)
-  (array-slice-ref (ensure-matrix 'submatrix a) (list row-range col-range)))
-} 
- TODO
+                 (define M (matrix ([1 2 3] [4 5 6])))
+                 (matrix-ref M 0 2)
+                 (matrix-ref M 1 2)]
 }
 
 @deftogether[(@defproc[(matrix-row [M (Matrix A)] [i Integer]) (Matrix A)]
               @defproc[(matrix-col [M (Matrix A)] [j Integer]) (Matrix A)])]{
-Returns the given row or column.
+Returns the @racket[i]th row or @racket[j]th column as a matrix.
 @examples[#:eval untyped-eval
-                 (define A (matrix ([1 2 3] [4 5 6])))
-                 (matrix-row A 1)
-                 (matrix-col A 0)]
+                 (define M (matrix ([1 2 3] [4 5 6])))
+                 (matrix-row M 1)
+                 (matrix-col M 0)]
+}
+
+@defproc[(submatrix [M (Matrix A)]
+                    [is (U Slice (Sequenceof Integer))]
+                    [js (U Slice (Sequenceof Integer))])
+         (Array A)]{
+Returns a submatrix or subarray of @racket[M], where @racket[is] and @racket[js] specify
+respectively the rows and columns to keep. Like @racket[array-slice-ref], but constrained so the
+result has exactly two axes.
+@examples[#:eval typed-eval
+                 (submatrix (identity-matrix 5) (:: 1 #f 2) (::))
+                 (submatrix (identity-matrix 5) '() '(1 2 4))]
+Note that @racket[submatrix] may return an empty array, which is not a @tech{matrix}.
 }
 
 @defproc[(matrix-diagonal [M (Matrix A)]) (Array A)]{
-Returns array of the elements on the diagonal of the square matrix.
+Returns array of the entries on the diagonal of @racket[M].
 @examples[#:eval untyped-eval
-                 (define A (matrix ([1 2 3] [4 5 6] [7 8 9])))
-                 (matrix-diagonal A)]        
+                 (matrix-diagonal
+                  (matrix ([1 2 3] [4 5 6] [7 8 9])))]
 }
 
-@deftogether[(@defproc[(matrix-upper-triangle [M (Matrix A)]) (Matrix A)]
-              @defproc[(matrix-lower-triangle [M (Matrix A)]) (Matrix A)])]{                                                           
+@deftogether[(@defproc[(matrix-upper-triangle [M (Matrix A)]) (Matrix (U A 0))]
+              @defproc[(matrix-lower-triangle [M (Matrix A)]) (Matrix (U A 0))])]{
 The function @racket[matrix-upper-triangle] returns an upper
 triangular matrix (entries below the diagonal are zero) with
-elements from the given matrix. Likewise the function 
-@racket[matrix-lower-triangle] returns an lower triangular
+entries from the given matrix. Likewise the function 
+@racket[matrix-lower-triangle] returns a lower triangular
 matrix.
 @examples[#:eval untyped-eval
-                 (define A (matrix ([1 2 3] [4 5 6] [7 8 9])))
-                 (matrix-upper-triangle A)
-                 (matrix-lower-triangle A)]
+                 (define M (array+ (array 1) (axis-index-array #(5 7) 1)))
+                 M
+                 (matrix-upper-triangle M)
+                 (matrix-lower-triangle M)]
 }
 
 @deftogether[(@defproc[(matrix-rows [M (Matrix A)]) (Listof (Matrix A))]
@@ -392,45 +478,61 @@ matrix.
 The functions respectively returns a list of the rows or columns
 of the matrix.
 @examples[#:eval untyped-eval
-                 (define A (matrix ([1 2 3] [4 5 6])))
-                 (matrix-rows A)
-                 (matrix-cols A)]
+                 (define M (matrix ([1 2 3] [4 5 6])))
+                 (matrix-rows M)
+                 (matrix-cols M)]
 }
 
 @deftogether[(@defproc[(matrix-augment [Ms (Listof (Matrix A))]) (Matrix A)]
               @defproc[(matrix-stack [Ms (Listof (Matrix A))]) (Matrix A)])]{
 The function @racket[matrix-augment] returns a matrix whose columns are
-the columns of the matrices in @racket[Ms]. This implies that the matrices 
+the columns of the matrices in @racket[Ms]. The matrices 
 in list must have the same number of rows.
 
 The function @racket[matrix-stack] returns a matrix whose rows are
-the rows of the matrices in @racket[Ms]. This implies that the matrices 
+the rows of the matrices in @racket[Ms]. The matrices 
 in list must have the same number of columns.
 @examples[#:eval untyped-eval
-                 (define A (matrix ([1 1] [1 1])))
-                 (define B (matrix ([2 2] [2 2])))
-                 (define C (matrix ([3 3] [3 3])))
-                 (matrix-augment (list A B C))
-                 (matrix-stack (list A B C))]
+                 (define M0 (matrix ([1 1] [1 1])))
+                 (define M1 (matrix ([2 2] [2 2])))
+                 (define M2 (matrix ([3 3] [3 3])))
+                 (matrix-augment (list M0 M1 M2))
+                 (matrix-stack (list M0 M1 M2))]
 }
 
 @deftogether[
-(@defproc[(matrix-map-rows 
-           [f ((Matrix A) -> (U #f (Matrix B)))] [M (Matrix A)] [fail (-> F) (λ () #f)]) (Matrix B)]  
- @defproc[(matrix-map-cols 
-           [f ((Matrix A) -> (U #f (Matrix B)))] [M (Matrix A)] [fail (-> F) (λ () #f)]) (Matrix B)])]{
-In the simple case the function @racket[matrix-map-rows] applies the function @racket[f]
-to each row of @racket[M]. If the rows are called @racket[r0], @racket[r1], ... then
-the result matrix has the rows @racket[(f r0)], @racket[(f r1)], ... .
-In the three argument case, the result of @racket[(fail)] is used, 
-if @racket[f] returns @racket[#f].
-
-The function @racket[matrix-map-cols] works likewise but on rows.
+(@defproc*[([(matrix-map-rows 
+              [f ((Matrix A) -> (Matrix B))] [M (Matrix A)]) (Matrix B)]
+            [(matrix-map-rows
+              [f ((Matrix A) -> (U #f (Matrix B)))] [M (Matrix A)] [fail (-> F)])
+             (U F (Matrix B))])])]{
+The two-argument case applies the function @racket[f] to each row of @racket[M].
+If the rows are called @racket[r0], @racket[r1], ..., the result matrix has the rows
+@racket[(f r0)], @racket[(f r1)], ....
 
 @examples[#:eval untyped-eval
-                 (define A (matrix ([1 2 3] [4 5 6] [7 8 9] [10 11 12])))
+                 (define M (matrix ([1 2 3] [4 5 6] [7 8 9] [10 11 12])))
                  (define (double-row r) (matrix-scale r 2))
-                 (matrix-map-rows double-row A)]
+                 (matrix-map-rows double-row M)]
+
+In the three argument case, if @racket[f] returns @racket[#f], the result of @racket[(fail)] is
+returned:
+@interaction[#:eval typed-eval
+                    (define Z (make-matrix 4 4 0))
+                    Z
+                    (matrix-map-rows (λ: ([r : (Matrix Real)])
+                                       (matrix-normalize r 2 (λ () #f)))
+                                     Z
+                                     (λ () 'FAILURE))]
+}
+
+@deftogether[
+(@defproc*[([(matrix-map-cols 
+              [f ((Matrix A) -> (Matrix B))] [M (Matrix A)]) (Matrix B)]
+            [(matrix-map-cols
+              [f ((Matrix A) -> (U #f (Matrix B)))] [M (Matrix A)] [fail (-> F)])
+             (U F (Matrix B))])])]{
+Like @racket[matrix-map-cols], but maps @racket[f] over columns.
 }
 
 
@@ -439,16 +541,16 @@ The function @racket[matrix-map-cols] works likewise but on rows.
 
 @section[#:tag "matrix:basic"]{Basic Operations}
 
-@defproc[(matrix-conjugate [M (Matrix A)]) (Matrix A)]{
-Returns a matrix where each element of the given matrix is conjugated.
+@defproc[(matrix-conjugate [M (Matrix Number)]) (Matrix Number)]{
+Returns a matrix where each entry of the given matrix is conjugated.
 @examples[#:eval untyped-eval
                  (matrix-conjugate (matrix ([1 +i] [-1 2+i])))]
 }
 
-@margin-note{@hyperlink["http://en.wikipedia.org/wiki/Transpose"]{Wikipedia: Transpose}}
+@margin-note*{@hyperlink["http://en.wikipedia.org/wiki/Transpose"]{Wikipedia: Transpose}}
 @deftogether[(@defproc[(matrix-transpose [M (Matrix A)]) (Matrix A)]
-              @defproc[(matrix-hermitian [M (Matrix A)]) (Matrix A)])]{
-@margin-note{@hyperlink["http://en.wikipedia.org/wiki/Hermitian_matrix"]{Wikipedia: Hermitian}}
+              @defproc[(matrix-hermitian [M (Matrix Number)]) (Matrix Number)])]{
+@margin-note*{@hyperlink["http://en.wikipedia.org/wiki/Hermitian_matrix"]{Wikipedia: Hermitian}}
 Returns the transpose or the hermitian of the matrix. 
 The hermitian of a matrix is the conjugate of the transposed matrix.
 For a real matrix these operations return the the same result.
@@ -457,71 +559,68 @@ For a real matrix these operations return the the same result.
                  (matrix-hermitian (matrix ([1 +i] [2 +2i] [3 +3i])))]
 }
 
-@margin-note{@hyperlink["http://en.wikipedia.org/wiki/Trace_(linear_algebra)"]{Wikipedia: Trace}}
-@defproc[(matrix-trace [M (Matrix Number)]) (Matrix Number)]{
+@margin-note*{@hyperlink["http://en.wikipedia.org/wiki/Trace_(linear_algebra)"]{Wikipedia: Trace}}
+@defproc[(matrix-trace [M (Matrix Number)]) Number]{
 Returns the trace of the square matrix. The trace of matrix is the 
-the sum of the diagonal elements. 
+the sum of the diagonal entries. 
 @examples[#:eval untyped-eval
                  (matrix-trace (matrix ([1 2] [3 4])))]
 }
 
 
-
 @;{==================================================================================================}
 
 
 @section[#:tag "matrix:inner"]{Inner Product Space Operations}
 
-The set of matrices of a given size forms a vector space.
-Since the vector space of @racket[mx1] matrices is isomorphic to
-the vector space of vectors of size @racket[m], any inner
-product (or norm) induce an inner product (or norm) on vectors.
-
-Put differently, the following innner products and norm, even
-though defined on general matrices, work on vectors in the
-form of column and row matrices.
+The following functions treat matrices as vectors in an inner product space.
+It often makes most sense to use these vector-space functions only for
+@tech[#:key "row matrix"]{row matrices} and @tech[#:key "column matrix"]{column matrices}, which are
+essentially vectors as we normally think of them.
+There are exceptions, however, such as the fact that the Frobenius or Euclidean norm (implemented by
+@racket[matrix-2norm]) can be used to measure error between matrices in a way that meets certain
+reasonable criteria (specifically, it is submultiplicative).
 
 See @secref{matrix:op-norm} for similar functions (e.g. norms and angles) defined by considering
 matrices as operators between inner product spaces consisting of column matrices.
 
-@margin-note{@hyperlink["http://en.wikipedia.org/wiki/Norm_(mathematics)"]{Wikipedia: Norm}}
-@deftogether[(@defproc[(matrix-1norm [M (Matrix Number)]) Number]
-              @defproc[(matrix-2norm [M (Matrix Number)]) Number]
-              @defproc[(matrix-inf-norm [M (Matrix Number)]) Number]
-              @defproc*[([(matrix-norm [M (Matrix Number)]) Number]
-                         [(matrix-norm [M (Matrix Number)] [p Real]) Number])])]{
-The first three functions compute the L1-norm, the L2-norm, and, the L∞-norm respectively.
+@margin-note*{@hyperlink["http://en.wikipedia.org/wiki/Norm_(mathematics)"]{Wikipedia: Norm}}
+@deftogether[(@defproc[(matrix-1norm [M (Matrix Number)]) Nonnegative-Real]
+              @defproc[(matrix-2norm [M (Matrix Number)]) Nonnegative-Real]
+              @defproc[(matrix-inf-norm [M (Matrix Number)]) Nonnegative-Real]
+              @defproc[(matrix-norm [M (Matrix Number)] [p Real 2]) Nonnegative-Real])]{
+Respectively compute the L@subscript{1} norm, L@subscript{2} norm, L@subscript{∞}, and
+L@subscript{p} norm.
 
-The L1-norm is also known under the names Manhattan- or taxicab-norm.
-The L1-norm of a matrix is the sum of magnitudes of the entries in the matrix.
+The L@subscript{1} norm is also known under the names Manhattan or taxicab norm.
+The L@subscript{1} norm of a matrix is the sum of magnitudes of the entries in the matrix.
 
-The L2-norm is also known under the names Euclidean- or Frobenius-norm. 
-The L2-norm of a matrix is the square root of the sum of squares of 
+The L@subscript{2} norm is also known under the names Euclidean or Frobenius norm.
+The L@subscript{2} norm of a matrix is the square root of the sum of squares of 
 magnitudes of the entries in the matrix.
 
-The L∞-norm is also known as the maximum- or infinity-norm.
-The L∞-norm computes the maximum magnitude of the entries in the matrix.
+The L@subscript{∞} norm is also known as the maximum or infinity norm.
+The L@subscript{∞} norm computes the maximum magnitude of the entries in the matrix.
 
-The function @racket[matrix-norm] computes the Lp-norm. 
-For a number @racket[p>=1] the @racket[p]th root of the sum
-of all entries to the @racket[p]th power. 
+For @racket[p >= 1], @racket[matrix-norm] computes the L@subscript{@racket[p]} norm: 
+the @racket[p]th root of the sum of all entry magnitudes to the @racket[p]th power.
 @;{MathJax would be nice to have in Scribble...}
-If no @racket[p] is given, the 2-norm (Eucledian) is used.
 @examples[#:eval untyped-eval
                  (matrix-1norm    (col-matrix [1 2]))
                  (matrix-2norm    (col-matrix [1 2]))
                  (matrix-inf-norm (col-matrix [1 2]))
-                 (matrix-norm     (col-matrix [1 2]) 3)]
+                 (matrix-norm     (col-matrix [1 2]) 3)
+                 (matrix-norm     (col-matrix [1 2]) +inf.0)]
 }
 
 @defproc*[([(matrix-dot [M (Matrix Number)]) Nonnegative-Real]
-           [(matrix-dot [M1 (Matrix Number)] [M2 (Matrix Number)]) Number])]{
+           [(matrix-dot [M (Matrix Number)] [N (Matrix Number)]) Number])]{
 
-The call @racket[(matrix-dot M1 M2)] computes the Frobenius inner product of the 
+The call @racket[(matrix-dot M N)] computes the Frobenius inner product of the 
 two matrices with the same shape.
 In other words the sum of @racket[(* a (conjugate b))] is computed where 
-@racket[a] runs over the entries in @racket[M1] and @racket[b] runs over
-the corresponding entries in @racket[M2].
+@racket[a] runs over the entries in @racket[M] and @racket[b] runs over
+the corresponding entries in @racket[N].
 
 The call @racket[(matrix-dot M)] computes @racket[(matrix-dot M M)] efficiently.
 @examples[#:eval untyped-eval
@@ -529,79 +628,130 @@ The call @racket[(matrix-dot M)] computes @racket[(matrix-dot M M)] efficiently.
                  (+ (* 1 3) (* 2 4))]
 }
 
-@defproc[(matrix-cos-angle [M (Matrix Number)]) Number]{
+@defproc[(matrix-cos-angle [M (Matrix Number)] [N (Matrix Number)]) Number]{
 Returns the cosine of the angle between two matrices w.r.t. the inner produce space induced by
 the Frobenius inner product. That is it returns 
 
 @racket[(/ (matrix-dot M N) (* (matrix-2norm M) (matrix-2norm N)))]
 
 @examples[#:eval untyped-eval
-                 (define e1 (col-matrix [1 0]))
-                 (define e2 (col-matrix [0 1]))
-                 (matrix-cos-angle e1 e2)
-                 (matrix-cos-angle e1 (matrix+ e1 e2))]
+                 (define M (col-matrix [1 0]))
+                 (define N (col-matrix [0 1]))
+                 (matrix-cos-angle M N)
+                 (matrix-cos-angle M (matrix+ M N))]
 }
 
-@defproc[(matrix-angle [M1 (Matrix Number)] [M2 (Matrix Number)]) Number]{
-Equivalent to @racket[(acos (matrix-cos-angle M0 M1))].
+@defproc[(matrix-angle [M (Matrix Number)] [N (Matrix Number)]) Number]{
+Equivalent to @racket[(acos (matrix-cos-angle M N))].
 @examples[#:eval untyped-eval
-                 (require racket/math) ; for radians->degrees
-                 (define e1 (col-matrix [1 0]))
-                 (define e2 (col-matrix [0 1]))
-                 (radians->degrees (matrix-angle e1 e2))
-                 (radians->degrees (matrix-angle e1 (matrix+ e1 e2)))]
+                 (require (only-in math/base radians->degrees))
+                 (define M (col-matrix [1 0]))
+                 (define N (col-matrix [0 1]))
+                 (radians->degrees (matrix-angle M N))
+                 (radians->degrees (matrix-angle M (matrix+ M N)))]
 }
 
-@defproc[(matrix-normalize [M (Matrix Number)] [p Real 2] [fail (-> A) raise-argument-error]) (U (Matrix Number) A)]{
-To normalize a matrix is to scale it, such that the result has norm 1.
-                                                                                                              
-The call @racket[(matrix-normalize M p fail)] normalizes @racket[M] with respect to 
-the @racket[Lp]-norm. If the matrix @racket[M] has norm 0, the result of calling 
-the thunk @racket[fail] is returned.
-
-If no fail-thunk is given, an argument-error exception is raised.
-If no @racket[p] the L2-norm (Euclidean norm) is used.
-@examples[#:eval untyped-eval
-                 (require racket/math) ; for radians->degrees
+@defproc[(matrix-normalize [M (Matrix Number)] [p Real 2] [fail (-> F) (λ () (error ...))])
+         (U F (Matrix Number))]{
+Normalizes @racket[M] with respect to the L@subscript{@racket[p]} norm.
+@examples[#:eval typed-eval
                  (matrix-normalize (col-matrix [1 1]))
                  (matrix-normalize (col-matrix [1 1]) 1)
                  (matrix-normalize (col-matrix [1 1]) +inf.0)]
+The result of applying the @tech{failure thunk} @racket[fail] is returned if @racket[M]'s norm is
+zero.
 }
 
-@deftogether[(@defproc[(matrix-normalize-rows [M (matrix Number)] [p Real 2] [fail (-> A) raise-argument-error]) (Matrix Number)]
-              @defproc[(matrix-normalize-cols [M (matrix Number)] [p Real 2] [fail (-> A) raise-argument-error]) (Matrix Number)])]{
+@deftogether[(@defproc[(matrix-normalize-rows [M (Matrix Number)]
+                                              [p Real 2]
+                                              [fail (-> F) (λ () (error ...))]) (Matrix Number)]
+              @defproc[(matrix-normalize-cols [M (Matrix Number)]
+                                              [p Real 2]
+                                              [fail (-> F) (λ () (error ...))]) (Matrix Number)])]{
 As @racket[matrix-normalize] but each row or column is normalized separately.
-The result is this a matrix with unit vectors as rows or columns.
-@examples[#:eval untyped-eval
-                 (require racket/math) ; for radians->degrees
+The result is a matrix with unit vectors as rows or columns.
+@examples[#:eval typed-eval
                  (matrix-normalize-rows (matrix [[1 2] [2 4]]))
                  (matrix-normalize-cols (matrix [[1 2] [2 4]]))]
+The result of applying the @tech{failure thunk} @racket[fail] is returned if the norm of any row
+or column in @racket[M] is zero.
 }
-                                                                                                             
-@deftogether[(@defproc[(matrix-rows-orthogonal? [M (Matrix Number)] [eps Real (* 10 epsilon.0)]) Boolean]
-              @defproc[(matrix-cols-orthogonal? [M (Matrix Number)] [eps Real (* 10 epsilon.0)]) Boolean])]{
+
+@deftogether[(@defproc[(matrix-rows-orthogonal? [M (Matrix Number)]
+                                                [eps Real (* 10 epsilon.0)])
+                       Boolean]
+              @defproc[(matrix-cols-orthogonal? [M (Matrix Number)]
+                                                [eps Real (* 10 epsilon.0)])
+                       Boolean])]{
 Returns @racket[#t] if the rows or columns of @racket[M] 
 are very close of being orthogonal (by default a few epsilons). 
 @examples[#:eval untyped-eval
-                 (require racket/math) ; for radians->degrees
                  (matrix-rows-orthogonal? (matrix [[1 1] [-1 1]]))
                  (matrix-cols-orthogonal? (matrix [[1 1] [-1 1]]))]
 }
 
 
-
 @;{==================================================================================================}
 
 
 @section[#:tag "matrix:solve"]{Solving Systems of Equations}
 
-@defthing[matrix-solve Procedure]{}
+@defproc[(matrix-solve [M (Matrix Number)] [B (Matrix Number)] [fail (-> F) (λ () (error ...))])
+         (U F (Matrix Number))]{
+Returns the matrix @racket[X] for which @racket[(matrix* M X)] is @racket[B].
+@racket[M] and @racket[B] must have the same number of rows.
 
-@defthing[matrix-inverse Procedure]{}
+It is typical for @racket[B] (and thus @racket[X]) to be a @tech{column matrix}, but not required.
+If @racket[B] is not a column matrix, @racket[matrix-solve] solves for all the columns in @racket[B]
+simultaneously.
 
-@defthing[matrix-invertible? Procedure]{}
+@examples[#:eval typed-eval
+                 (define M (matrix [[7 5] [3 -2]]))
+                 (define B0 (col-matrix [3 22]))
+                 (define B1 (col-matrix [19 4]))
+                 (matrix-solve M B0)
+                 (matrix* M (col-matrix [4 -5]))
+                 (matrix-solve M B1)
+                 (matrix-cols (matrix-solve M (matrix-augment (list B0 B1))))]
 
-@defthing[matrix-determinant Procedure]{}
+@racket[matrix-solve] does not solve overconstrained or underconstrained systems, meaning that
+@racket[M] must be invertible.
+If @racket[M] is not invertible, the result of applying the @tech{failure thunk} @racket[fail] is
+returned.
+
+@racket[matrix-solve] is implemented using @racket[matrix-gauss-elim] to preserve exactness in its
+output, with partial pivoting for greater numerical stability when @racket[M] is not exact.
+
+See @racket[vandermonde-matrix] for an example that uses @racket[matrix-solve] to compute Legendre
+polynomials.
+}
+
+@define[inverse-url]{http://en.wikipedia.org/wiki/Invertible_matrix}
+@margin-note*{@hyperlink[inverse-url]{Wikipedia: Invertible Matrix}}
+@defproc[(matrix-inverse [M (Matrix Number)] [fail (-> F) (λ () (error ...))]) (U F (Matrix Number))]{
+Returns the inverse of @racket[M] if it exists; otherwise returns the result of applying the
+@tech{failure thunk} @racket[fail].
+@examples[#:eval typed-eval
+                 (matrix-inverse (identity-matrix 3))
+                 (matrix-inverse (matrix [[7 5] [3 -2]]))
+                 (matrix-inverse (matrix [[1 2] [10 20]]))
+                 (matrix-inverse (matrix [[1 2] [10 20]]) (λ () #f))]
+}
+
+@defproc[(matrix-invertible? [M (Matrix Number)]) Boolean]{
+Returns @racket[#t] when @racket[M] is a @racket[square-matrix?] and @racket[(matrix-determinant M)]
+is nonzero.
+}
+
+@margin-note*{@hyperlink["http://en.wikipedia.org/wiki/Determinant"]{Wikipedia: Determinant}}
+@defproc[(matrix-determinant [M (Matrix Number)]) Number]{
+Returns the determinant of @racket[M], which must be a @racket[square-matrix?].
+@examples[#:eval typed-eval
+                 (matrix-determinant (diagonal-matrix '(1 2 3 4)))
+                 (* 1 2 3 4)
+                 (matrix-determinant (matrix [[1 2] [10 20]]))
+                 (matrix-determinant (col-matrix [1 2]))]
+}
 
                                                                                   
 @;{==================================================================================================}
@@ -609,11 +759,87 @@ are very close of being orthogonal (by default a few epsilons).
 
 @section[#:tag "matrix:row-alg"]{Row-Based Algorithms}
 
-@defthing[matrix-gauss-elim Procedure]{}
+@define[gauss-url]{http://en.wikipedia.org/wiki/Gaussian_elimination}
+@define[gauss-jordan-url]{http://en.wikipedia.org/wiki/Gauss%E2%80%93Jordan_elimination}
 
-@defthing[matrix-row-echelon Procedure]{}
+@margin-note*{@hyperlink[gauss-url]{Wikipedia: Gaussian elimination}}
+@defproc[(matrix-gauss-elim [M (Matrix Number)]
+                            [jordan? Any #f]
+                            [unitize-pivot? Any #f]
+                            [pivoting (U 'first 'partial) 'partial])
+         (Values (Matrix Number) (Listof Index))]{
+@margin-note*{@hyperlink[gauss-jordan-url]{Wikipedia: Gauss-Jordan elimination}}
+Implements Gaussian elimination or Gauss-Jordan elimination.
 
-@defthing[matrix-lu Procedure]{}
+If @racket[jordan?] is true, row operations are done both above and below the pivot.
+If @racket[unitize-pivot?] is true, the pivot's row is scaled so that the pivot value is @racket[1].
+When both are true, the algorithm is called Gauss-Jordan elimination, and the result matrix is in
+@italic{reduced} row echelon form.
+
+If @racket[pivoting] is @racket['first], the first nonzero entry in the current column is used as the
+pivot. If @racket[pivoting] is @racket['partial], the largest-magnitude nonzero entry is used, which
+improves numerical stability on average when @racket[M] contains inexact entries.
+
+The first return value is the result of Gaussian elimination.
+
+The second return value is a list of indexes of columns that did not have a nonzero pivot.
+
+See @racket[matrix-row-echelon] for examples.
+}
+
+@define[row-echelon-url]{http://en.wikipedia.org/wiki/Row_echelon_form}
+
+@margin-note*{@hyperlink[row-echelon-url]{Wikipedia: Row echelon form}}
+@defproc[(matrix-row-echelon [M (Matrix Number)]
+                             [jordan? Any #f]
+                             [unitize-pivot? Any #f]
+                             [pivoting (U 'first 'partial) 'partial])
+         (Matrix Number)]{
+Like @racket[matrix-gauss-elim], but returns only the result of Gaussian elimination.
+@examples[#:eval typed-eval
+                 (define M (matrix [[2 1 -1] [-3 -1 2] [-2 1 2]]))
+                 (matrix-row-echelon M)
+                 (matrix-row-echelon M #t)
+                 (matrix-identity? (matrix-row-echelon M #t #t))]
+The last example shows that @racket[M] is invertible.
+
+Using @racket[matrix-row-echelon] to solve a system of linear equations (without checking for
+failure):
+@interaction[#:eval typed-eval
+                    (define B (col-matrix [8 -11 -3]))
+                    (define MB (matrix-augment (list M B)))
+                    (matrix-col (matrix-row-echelon MB #t #t) 3)
+                    (matrix-solve M B)]
+
+Using @racket[matrix-row-echelon] to invert a matrix (also without checking for failure):
+@interaction[#:eval typed-eval
+                    (define MI (matrix-augment (list M (identity-matrix 3))))
+                    (submatrix (matrix-row-echelon MI  #t #t) (::) (:: 3 #f))
+                    (matrix-inverse M)]
+}
+
+@define[lu-url]{http://en.wikipedia.org/wiki/LU_decomposition}
+@margin-note*{@hyperlink[lu-url]{Wikipedia: LU decomposition}}
+@defproc[(matrix-lu [M (Matrix Number)] [fail (-> F) (λ () (error ...))])
+         (Values (U F (Matrix Number)) (Matrix Number))]{
+Returns the LU decomposition of @racket[M] (which must be a @racket[square-matrix?]) if one exists.
+An LU decomposition exists if @racket[M] can be put in row-echelon form without swapping rows.
+
+Because @racket[matrix-lu] returns a @italic{unit} lower-triangular matrix (i.e. a lower-triangular
+matrix with only ones on the diagonal), the decomposition is unique if it exists.
+
+@examples[#:eval typed-eval
+                 (define-values (L U)
+                   (matrix-lu (matrix [[4 3] [6 3]])))
+                 (values L U)
+                 (matrix* L U)]
+
+If @racket[M] does not have an LU decomposition, the first result is the result of applying the
+@tech{failure thunk} @racket[fail], and the second result is the original argument @racket[M]:
+@interaction[#:eval typed-eval
+                    (matrix-lu (matrix [[0 1] [1 1]]))
+                    (matrix-lu (matrix [[0 1] [1 1]]) (λ () #f))]
+}
 
                                                                                   
 @;{==================================================================================================}
@@ -621,57 +847,143 @@ are very close of being orthogonal (by default a few epsilons).
 
 @section[#:tag "matrix:ortho-alg"]{Orthogonal Algorithms}
 
-@defthing[matrix-gram-schmidt Procedure]{}
+@define[gram-schmidt-url]{http://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process}
+@define[reortho-pdf]{http://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf}
 
-@defthing[matrix-basis-extension Procedure]{}
+@margin-note*{@hyperlink[gram-schmidt-url]{Wikipedia: Gram-Schmidt process}}
+@defproc[(matrix-gram-schmidt [M (Matrix Number)] [normalize? Any #f] [start-col Integer 0])
+         (Array Number)]{
+Returns an array whose columns are orthogonal and span the same subspace as @racket[M]'s columns.
+The number of columns in the result is the rank of @racket[M].
+If @racket[normalize?] is true, the columns are also normalized.
 
-@margin-note{@hyperlink["http://en.wikipedia.org/wiki/QR_decomposition"]{Wikipedia: QR decomposition}}
-@defproc*[([(matrix-qr [M (Matrix Real)])   (Values (Matrix Real) (Matrix Real))]
-           [(matrix-qr [M (Matrix Real)] [full Any]) (Values (Matrix Real) (Matrix Real))] 
-           [(matrix-qr [M (Matrix Number)])  (Values (Matrix Number) (Matrix Number))]
-           [(matrix-qr [M (Matrix Number)] [full Any]) (Values (Matrix Number) (Matrix Number))])]{
+@examples[#:eval typed-eval
+                 (define M 
+                   (matrix [[12 -51   4]
+                            [ 6 167 -68]
+                            [-4  24 -41]]))
+                 (matrix-gram-schmidt M)
+                 (matrix-gram-schmidt M #t)
+                 (matrix-cols-orthogonal? (matrix-gram-schmidt M))
+                 (matrix-orthonormal? (matrix-gram-schmidt M #t))]
+
+Examples with rank-deficient matrices:
+@interaction[#:eval typed-eval
+                    (matrix-gram-schmidt (matrix [[ 1 -2 1]
+                                                  [-2  4 9]
+                                                  [ 3 -6 7]]))
+                    (matrix-gram-schmidt (make-matrix 3 3 0))]
+
+When @racket[start-col] is positive, the Gram-Schmidt process is begun on column @racket[start-col]
+(but still using the previous columns to orthogonalize the remaining columns).
+This feature is generally not directly useful, but is used in the implementation of
+@racket[matrix-basis-extension].
+
+@margin-note*{* @bold{On the round-off error analysis of the Gram-Schmidt algorithm with
+                      reorthogonalization.},
+              Luc Giraud, Julien Langou and Miroslav Rozloznik. 2002. @hyperlink[reortho-pdf]{(PDF)}}
+While Gram-Schmidt with inexact matrices is known to be unstable, using it twice tends to remove
+instabilities:*
+@interaction[#:eval typed-eval
+                    (define M (matrix [[0.70000 0.70711]
+                                       [0.70001 0.70711]]))
+                    (matrix-orthonormal?
+                     (matrix-gram-schmidt M #t))
+                    (matrix-orthonormal?
+                     (matrix-gram-schmidt (matrix-gram-schmidt M) #t))]
+}
+
+@defproc[(matrix-basis-extension [M (Matrix Number)]) (Array Number)]{
+Returns additional orthogonal columns which, if augmented with @racket[M], would result in an
+orthogonal matrix of full rank. If @racket[M]'s columns are normalized, the result's columns are
+normalized.
+}
+
+@define[qr-url]{http://en.wikipedia.org/wiki/QR_decomposition}
+
+@margin-note*{@hyperlink["http://en.wikipedia.org/wiki/QR_decomposition"]{Wikipedia: QR decomposition}}
+@defproc*[([(matrix-qr [M (Matrix Number)])  (Values (Matrix Number) (Matrix Number))]
+           [(matrix-qr [M (Matrix Number)] [full? Any]) (Values (Matrix Number) (Matrix Number))])]{
 Computes a QR-decomposition of the matrix @racket[M]. The values returned are
-the matrices @racket[Q] and @racket[R]. If @racket[full] is false, then
+the matrices @racket[Q] and @racket[R]. If @racket[full?] is @racket[#f], then
 a reduced decomposition is returned, otherwise a full decomposition is returned.
 
-@margin-note{An @italic{orthonormal} matrix has columns which are orthooginal, unit vectors.}
+@margin-note*{An @italic{orthonormal} matrix has columns which are orthoginal, unit vectors.}
 The (full) decomposition of a square matrix consists of two matrices: 
 a orthogonal matrix @racket[Q] and an upper triangular matrix @racket[R], 
 such that @racket[QR = M]. 
 
 For tall non-square matrices @racket[R], the triangular part of the full decomposition,
 contains zeros below the diagonal. The reduced decomposition leaves the zeros out.
-See the Wikipedia entry on @hyperlink["http://en.wikipedia.org/wiki/QR_decomposition"]{QR decomposition}
+See the Wikipedia entry on @hyperlink[qr-url]{QR decomposition}
 for more details.
 
+@examples[#:eval typed-eval
+                 (define M
+                   (matrix [[12 -51   4]
+                            [ 6 167 -68]
+                            [-4  24 -41]]))
+                 (define-values (Q R) (matrix-qr M))
+                 (values Q R)
+                 (matrix= M (matrix* Q R))]
+
+The difference between full and reduced decompositions:
+@interaction[#:eval typed-eval
+                    (define M
+                      (matrix [[12 -51]
+                               [ 6 167]
+                               [-4  24]]))
+                    (define-values (Q1 R1) (matrix-qr M #f))
+                    (define-values (Q2 R2) (matrix-qr M #t))
+                    (values Q1 R1)
+                    (values Q2 R2)
+                    (matrix= M (matrix* Q1 R1))
+                    (matrix= M (matrix* Q2 R2))]
+
 The decomposition @racket[M = QR] is useful for solving the equation @racket[Mx=v]. 
 Since the inverse of Q is simply the transpose of Q, 
   @racket[Mx=v  <=>  QRx=v  <=>  Rx = Q^T v].
 And since @racket[R] is upper triangular, the system can be solved by back substitution.
 
 The algorithm used is Gram-Schmidt with reorthogonalization.
-See the paper @hyperlink["http://www.cerfacs.fr/algor/reports/2002/TR_PA_02_33.pdf"]{On the round-off error analysis of the Gram-Schmidt algorithm with reorthogonalization.}
-by Luc Giraud, Julien Langou, and, Miroslav Rozloznik.
-}                                                                                  
+}
+
+
 @;{==================================================================================================}
 
 
+@define[op-norm-url]{http://en.wikipedia.org/wiki/Matrix_norm#Induced_norm}
+
 @section[#:tag "matrix:op-norm"]{Operator Norms and Comparing Matrices}
 
+@secref{matrix:inner} describes functions that deal with matrices as vectors in
+an inner product space. This section describes functions that deal with matrices as
+@italic{linear operators}, or as functions from column matrices to column matrices.
+
+@margin-note*{@hyperlink[op-norm-url]{Wikipedia: Induced norm}}
+In this setting, a norm is the largest relative change in magnitude an operator (i.e. matrix) can
+effect on a column matrix, where ``magnitude'' is defined by a vector norm.
+(See the Wikipedia article linked to in the margin for a formal definition.)
+Matrix norms that are defined in terms of a vector norm are called @italic{induced norms}, or
+@deftech{operator norms}.
+
 @defproc[(matrix-op-1norm [M (Matrix Number)]) Nonnegative-Real]{
-TODO: describe
+The @tech{operator norm} induced by the vector norm @racket[matrix-1norm].
 
 When M is a column matrix, @racket[(matrix-op-1norm M)] is equivalent to @racket[(matrix-1norm M)].
 }
 
 @defproc[(matrix-op-2norm [M (Matrix Number)]) Nonnegative-Real]{
-TODO: describe (spectral norm)
+The @tech{operator norm} induced by the vector norm @racket[matrix-2norm].
+                                             
+This function is currently undefined because a required algorithm (singular value decomposition or
+eigendecomposition) is not yet implemented in @racketmodname[math/matrix].
 
 When M is a column matrix, @racket[(matrix-op-2norm M)] is equivalent to @racket[(matrix-2norm M)].
 }
 
 @defproc[(matrix-op-inf-norm [M (Matrix Number)]) Nonnegative-Real]{
-TODO: describe
+The @tech{operator norm} induced by the vector norm @racket[matrix-inf-norm].
 
 When M is a column matrix, @racket[(matrix-op-inf-norm M)] is equivalent to
 @racket[(matrix-inf-norm M)].
@@ -680,16 +992,22 @@ When M is a column matrix, @racket[(matrix-op-inf-norm M)] is equivalent to
 @defproc[(matrix-basis-cos-angle [M0 (Matrix Number)] [M1 (Matrix Number)]) Number]{
 Returns the cosine of the angle between the two subspaces spanned by @racket[M0] and @racket[M1].
 
+This function is currently undefined because a required algorithm (singular value decomposition or
+eigendecomposition) is not yet implemented in @racketmodname[math/matrix].
+
 When @racket[M0] and @racket[M1] are column matrices, @racket[(matrix-basis-cos-angle M0 M1)] is
 equivalent to @racket[(matrix-cos-angle M0 M1)].
 }
          
 @defproc[(matrix-basis-angle [M0 (Matrix Number)] [M1 (Matrix Number)]) Number]{
 Equivalent to @racket[(acos (matrix-basis-cos-angle M0 M1))].
+
+The function is currently undefined because @racket[matrix-basis-cos-angle] is currently undefined.
 }
 
 @defparam[matrix-error-norm norm ((Matrix Number) -> Nonnegative-Real)]{
 The norm used by @racket[matrix-relative-error] and @racket[matrix-absolute-error].
+The default value is @racket[matrix-op-inf-norm].
 
 Besides being a true norm, @racket[norm] should also be @deftech{submultiplicative}:
 @racketblock[(norm (matrix* M0 M1)) <= (* (norm A) (norm B))]
@@ -740,7 +1058,7 @@ Equivalent to
 @defproc[(matrix-orthonormal? [M (Matrix Number)] [eps Real (* 10 epsilon.0)]) Boolean]{
 Returns @racket[#t] when @racket[M] is very close to being orthonormal; that is, when
 @racket[(matrix* M (matrix-hermitian M))] is very close to an identity matrix.
-In fact, @racket[(matrix-orthonormal? M eps)] is equivalent to
+Equivalent to
 @racketblock[(matrix-identity? (matrix* M (matrix-hermitian M)) eps)]
 }
 
diff --git a/collects/math/scribblings/math.scrbl b/collects/math/scribblings/math.scrbl
index 9b00949c96..02338bddf4 100644
--- a/collects/math/scribblings/math.scrbl
+++ b/collects/math/scribblings/math.scrbl
@@ -17,17 +17,18 @@ for working with numbers and collections of numbers. These include
   @item{@racketmodname[math/bigfloat]: Arbitrary-precision floating-point functions}
   @item{@racketmodname[math/number-theory]: Number-theoretic functions}
   @item{@racketmodname[math/array]: Functional arrays for operating on large rectangular data sets}
-  @item{@racketmodname[math/matrix]: Linear algebra functions for arrays @bold{(currently being documented)}}
+  @item{@racketmodname[math/matrix]: Linear algebra functions for arrays}
   @item{@racketmodname[math/distributions]: Probability distributions}
   @item{@racketmodname[math/statistics]: Statistical functions}
 ]
 
 With this library, we hope to support a wide variety of applied mathematics in
 Racket, including simulation, statistical inference, signal processing, and
-combinatorics. If you find it lacking for your variety of mathematics,
-please contact us or post to one of the
-@hyperlink["http://racket-lang.org/community.html"]{mailing lists}
-to make suggestions or submit patches.
+combinatorics. If you find it lacking for your variety of mathematics, please
+@itemlist[@item{Visit the @hyperlink["https://github.com/plt/racket/wiki/Math-Library-Features"]{Math Library Features wiki page}
+                to see what is planned.}
+          @item{Contact us or post to one of the @hyperlink["http://racket-lang.org/community.html"]{mailing lists} to
+                make suggestions or submit patches.}]
 
 @bold{This is a Typed Racket library.} It is most efficient to use it in Typed
 Racket, so that contracts are checked statically. However, almost all of it can
diff --git a/collects/math/tests/matrix-strictness-tests.rkt b/collects/math/tests/matrix-strictness-tests.rkt
index 5b2b61f264..697197d48b 100644
--- a/collects/math/tests/matrix-strictness-tests.rkt
+++ b/collects/math/tests/matrix-strictness-tests.rkt
@@ -96,7 +96,7 @@
     (check-true (array-strict? (matrix-col M i))))
   
   (for*: ([M  (list nonstrict-2x2-arr strict-2x2-arr)]
-          [spec  (list '(0) 0)])
+          [spec  (list '(0) (:: #f #f 2))])
     (check-true (array-strict? (submatrix M (::) spec))))
   )
 
@@ -153,7 +153,7 @@
       (check-false (array-strict? R))))
   
   (for*: ([M  (list nonstrict-2x2-arr strict-2x2-arr)]
-          [spec  (list '(0) 0)])
+          [spec  (list '(0) (:: #f #f 2))])
     (check-false (array-strict? (submatrix M (::) spec))))
   
   (for*: ([M  (list nonstrict-2x2-arr strict-2x2-arr)]