diff --git a/sparse-linear-algebra.cabal b/sparse-linear-algebra.cabal
--- a/sparse-linear-algebra.cabal
+++ b/sparse-linear-algebra.cabal
@@ -1,6 +1,6 @@
 name:                sparse-linear-algebra
-version:             0.2.9.1
-synopsis:            Sparse linear algebra in native Haskell.
+version:             0.2.9.2
+synopsis:            Numerical computation in native Haskell
 description:
   /Overview/
   .
diff --git a/src/Numeric/LinearAlgebra/Class.hs b/src/Numeric/LinearAlgebra/Class.hs
--- a/src/Numeric/LinearAlgebra/Class.hs
+++ b/src/Numeric/LinearAlgebra/Class.hs
@@ -37,11 +37,13 @@
 
 
 
--- * Matrix and vector elements (possibly Complex)
+-- * Matrix and vector elements (optionally Complex)
 class (Eq e , Fractional e, Floating e, Num (EltMag e), Ord (EltMag e)) => Elt e where
   type EltMag e :: *
+  -- | Complex conjugate, or identity function if its input is real-valued
   conj :: e -> e
   conj = id
+  -- | Magnitude
   mag :: e -> EltMag e
 
 instance Elt Double where {type EltMag Double = Double ; mag = id}
@@ -68,8 +70,8 @@
 v ./ n = recip n .* v
 
 -- | Convex combination of two vectors (NB: 0 <= `a` <= 1). 
-lerp :: (VectorSpace e, Num (Scalar e)) => Scalar e -> e -> e -> e
-lerp a u v = a .* u ^+^ ((1-a) .* v)
+cvx :: (VectorSpace e, Num (Scalar e)) => Scalar e -> e -> e -> e
+cvx a u v = a .* u ^+^ ((1-a) .* v)
 
 
 -- linearCombination :: (VectorSpace v , Foldable t) => t (Scalar v, v) -> v
@@ -110,18 +112,27 @@
 class (InnerSpace v, Num (RealScalar v), Eq (RealScalar v), Epsilon (Magnitude v), Show (Magnitude v), Ord (Magnitude v)) => Normed v where
   type Magnitude v :: *
   type RealScalar v :: *
-  norm1 :: v -> Magnitude v      -- ^ L1 norm
-  norm2Sq :: v -> Magnitude v    -- ^ Euclidean (L2) norm
-  normP :: RealScalar v -> v -> Magnitude v -- ^ Lp norm (p > 0)
-  normalize :: RealScalar v -> v -> v  -- ^ Normalize w.r.t. Lp norm
-  normalize2 :: v -> v     -- ^ Normalize w.r.t. L2 norm
-  normalize2' :: Floating (Scalar v) => v -> v -- ^ Normalize w.r.t. norm2' instead of norm2
+  -- | L1 norm
+  norm1 :: v -> Magnitude v
+  -- | Euclidean (L2) norm squared
+  norm2Sq :: v -> Magnitude v
+  -- | Lp norm (p > 0)
+  normP :: RealScalar v -> v -> Magnitude v
+  -- | Normalize w.r.t. Lp norm
+  normalize :: RealScalar v -> v -> v
+  -- | Normalize w.r.t. L2 norm
+  normalize2 :: v -> v
+  -- | Normalize w.r.t. norm2' instead of norm2
+  normalize2' :: Floating (Scalar v) => v -> v 
   normalize2' x = x ./ norm2' x
-  norm2 :: Floating (Magnitude v) => v -> Magnitude v -- ^ Euclidean norm
+  -- | Euclidean (L2) norm
+  norm2 :: Floating (Magnitude v) => v -> Magnitude v 
   norm2 x = sqrt (norm2Sq x)
-  norm2' :: Floating (Scalar v) => v -> Scalar v -- ^ Euclidean norm; returns a Complex (norm :+ 0) for containers of complex values
+  -- | Euclidean (L2) norm; returns a Complex (norm :+ 0) for Complex-valued vectors
+  norm2' :: Floating (Scalar v) => v -> Scalar v 
   norm2' x = sqrt $ x <.> x
-  norm :: Floating (Magnitude v) => RealScalar v -> v -> Magnitude v -- ^ Lp norm (p > 0)
+  -- | Lp norm (p > 0)
+  norm :: Floating (Magnitude v) => RealScalar v -> v -> Magnitude v 
   norm p v
     | p == 1 = norm1 v
     | p == 2 = norm2 v
@@ -205,11 +216,16 @@
 
 class (AdditiveGroup m, Epsilon (MatrixNorm m)) => MatrixRing m where
   type MatrixNorm m :: *
+  -- | Matrix-matrix product
   (##) :: m -> m -> m
-  (##^) :: m -> m -> m   -- ^ A B^T
-  (#^#) :: m -> m -> m   -- ^ A^T B
+  -- | A B^T
+  (##^) :: m -> m -> m
+  -- | A^T B
+  (#^#) :: m -> m -> m
   a #^# b = transpose a ## b
+  -- | Matrix transpose
   transpose :: m -> m
+  -- | Frobenius norm
   normFrobenius :: m -> MatrixNorm m
 
 
@@ -230,7 +246,9 @@
 
 class (VectorSpace v, MatrixRing (MatrixType v)) => LinearVectorSpace v where
   type MatrixType v :: *
+  -- | Matrix-vector action
   (#>) :: MatrixType v -> v -> v
+  -- | Dual matrix-vector action
   (<#) :: v -> MatrixType v -> v
 
 
@@ -247,7 +265,11 @@
 -- ** Linear systems
   
 class LinearVectorSpace v => LinearSystem v where
-  (<\>) :: (MonadIO m, MonadThrow m) => MatrixType v -> v -> m v
+  -- | Solve a linear system
+  (<\>) :: (MonadIO m, MonadThrow m) =>
+           MatrixType v   -- ^ System matrix
+        -> v              -- ^ Right-hand side
+        -> m v            -- ^ Result
 
 
 
@@ -262,6 +284,7 @@
 
 class Functor f => FiniteDim f where
   type FDSize f :: *
+  -- | Dimension (i.e. Int for SpVector, (Int, Int) for SpMatrix)
   dim :: f a -> FDSize f
 
 class FiniteDim' f where
@@ -301,6 +324,7 @@
 
 class HasData f a where
   type HDData f a :: *
+  -- | Number of nonzeros
   nnz :: f a -> Int
   dat :: f a -> HDData f a
 
@@ -313,6 +337,7 @@
 -- * Sparse : sparse datastructures
 
 class (FiniteDim f, HasData f a) => Sparse f a where
+  -- | Sparsity (fraction of nonzero elements)
   spy :: Fractional b => f a -> b
 
 class (FiniteDim' f, HasData' f) => Sparse' f where
@@ -323,10 +348,10 @@
 -- * Set : types that behave as sets
 
 class Functor f => Set f where
-  -- |union binary lift : apply function on _union_ of two "sets"
+  -- | Union binary lift : apply function on _union_ of two "sets"
   liftU2 :: (a -> a -> a) -> f a -> f a -> f a
 
-  -- |intersection binary lift : apply function on _intersection_ of two "sets"
+  -- | Intersection binary lift : apply function on _intersection_ of two "sets"
   liftI2 :: (a -> a -> b) -> f a -> f a -> f b
 
 
diff --git a/src/Numeric/LinearAlgebra/Sparse.hs b/src/Numeric/LinearAlgebra/Sparse.hs
--- a/src/Numeric/LinearAlgebra/Sparse.hs
+++ b/src/Numeric/LinearAlgebra/Sparse.hs
@@ -12,7 +12,7 @@
        (
          -- * Linear solvers
          -- ** Iterative methods
-         linSolve0, LinSolveMethod(..), (<\>),
+         (<\>),
          -- ** Moore-Penrose pseudoinverse
          pinv,
          -- ** Preconditioners
@@ -25,7 +25,8 @@
          triUpperSolve,
          -- * Eigensolvers
          eigsQR,
-         eigRayleigh,
+         -- eigRayleigh,
+         eigsArnoldi,
          -- * Matrix factorization algorithms
          -- ** QR
          qr,
@@ -43,7 +44,7 @@
          -- ** Condition number
          conditionNumberSM,
          -- ** Householder reflection
-         hhMat, hhRefl,
+         hhRefl,
          -- -- * Householder bidiagonalization         
          -- -- * Random arrays
          -- randArray,
@@ -51,14 +52,36 @@
          -- randMat, randVec, 
          -- -- ** Sparse "
          -- randSpMat, randSpVec,
-         -- * From/to SpVector
+         -- * Creation and conversion of sparse data
+         -- ** From/to SpVector
          fromListSV, toListSV,
-         -- * From/to SpMatrix
+         -- ** From/to SpMatrix
          fromListSM, toListSM,
+         -- * Operators
+         -- ** Inner product
+         (<.>),
+         -- ** Matrix-vector products
+         (#>), (<#),
+         -- ** Matrix-matrix products
+         (##), (#^#), (##^),
+         -- *** Sparsifying matrix-matrix products
+         (#~#), (#~^#), (#~#^),
+         -- ** Vector outer product
+         (><),
+         -- * Common operations
+         dim, nnz, spy,
+         -- ** Vector-related
+         (.*), (./), cvx,
+         norm, norm2, norm2', normalize, normalize2, normalize2',
+         norm1, hilbertDistSq,
+         -- ** Matrix-related
+         transpose, normFrobenius, 
          -- * Iteration combinators
          untilConvergedG0, untilConvergedG, untilConvergedGM,
          modifyInspectGuarded, modifyInspectGuardedM, IterationConfig (..),
-         modifyUntil, modifyUntilM
+         modifyUntil, modifyUntilM,
+         -- * Internal
+         linSolve0, LinSolveMethod(..)
        )
        where
 
@@ -101,7 +124,7 @@
 
 
 -- | A lumped constraint for numerical types
-type Num' x = (Epsilon x, Elt x, Show x, Ord x)
+type Num' x = (Epsilon x, Elt x, Show x, Ord x, Typeable x)
 
 
 
@@ -109,7 +132,7 @@
 
 -- * Matrix condition number
 
--- |uses the R matrix from the QR factorization
+-- | Matrix condition number: computes the QR factorization and extracts the extremal eigenvalues from the R factor
 conditionNumberSM :: (MonadThrow m, MatrixRing (SpMatrix a), Num' a, Typeable a) =>
      SpMatrix a -> m a
 conditionNumberSM m = do
@@ -135,10 +158,9 @@
   n = dim x
 
 
--- | Householder reflection: a vector `x` uniquely defines an orthogonal plane; the Householder operator reflects any point `v` with respect to this plane:
--- v' = (I - 2 x >< x) v
+-- | Householder reflection: a vector `x` uniquely defines an orthogonal (hyper)plane, i.e. an orthogonal subspace; the Householder operator reflects any point `v` through this subspace: v' = (I - 2 x >< x) v
 hhRefl :: Num a => SpVector a -> SpMatrix a
-hhRefl = hhMat (fromInteger 2)
+hhRefl = hhMat 2
 
 
 
@@ -247,7 +269,8 @@
 -- | Given a matrix A, returns a pair of matrices (Q, R) such that Q R = A, where Q is orthogonal and R is upper triangular. Applies Givens rotation iteratively to zero out sub-diagonal elements.
 {-# inline qr #-}
 qr :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, MonadThrow m) =>
-     SpMatrix a -> m (SpMatrix a, SpMatrix a)
+     SpMatrix a
+     -> m (SpMatrix a, SpMatrix a)  -- ^ Q, R
 qr mm = do 
      (qt, r, _) <- MTS.execStateT (modifyUntilM haltf qrstepf) gminit
      return (transpose qt, r) 
@@ -276,7 +299,7 @@
 -- ** QR algorithm
 
 -- | `eigsQR n mm` performs `n` iterations of the QR algorithm on matrix `mm`, and returns a SpVector containing all eigenvalues
-eigsQR :: (MonadThrow m, MonadIO m, Elt a, Normed (SpVector a), MatrixRing (SpMatrix a), Epsilon a, Typeable (Magnitude (SpVector a)), Typeable a, Show a) =>
+eigsQR :: (MonadThrow m, MonadIO m, Num' a, Normed (SpVector a), MatrixRing (SpMatrix a), Typeable (Magnitude (SpVector a))) =>
         Int
      -> Bool           -- ^ Print debug information        
      -> SpMatrix a     -- ^ Operand matrix
@@ -311,15 +334,14 @@
       return (b', mu')
 
 
--- | `eigArnoldi n aa b` computes at most n iterations of the Arnoldi algorithm to find a Krylov subspace of (A, b), along with a Hessenberg matrix of coefficients H. After that, it computes the QR decomposition of H, and the eigenvalues of A are listed on the diagonal of the R factor.
-eigArnoldi :: (Scalar (SpVector t) ~ t, MatrixType (SpVector t) ~ SpMatrix t,
-      Elt t, Normed (SpVector t), MatrixRing (SpMatrix t),
-      LinearVectorSpace (SpVector t), Epsilon t, MonadThrow m) =>
+-- | `eigsArnoldi n aa b` computes at most n iterations of the Arnoldi algorithm to find a Krylov subspace of (A, b), denoted Q, along with a Hessenberg matrix of coefficients H. After that, it computes the QR decomposition of H, denoted (O, R) and the eigenvalues of A are listed on the diagonal of the R factor.
+eigsArnoldi :: (Scalar (SpVector t) ~ t, MatrixType (SpVector t) ~ SpMatrix t,
+      Elt t, V (SpVector t), MatrixRing (SpMatrix t), Epsilon t, MonadThrow m) =>
      Int
      -> SpMatrix t
      -> SpVector t
      -> m (SpMatrix t, SpMatrix t, SpVector t) -- ^ Q, O, R
-eigArnoldi nitermax aa b = do
+eigsArnoldi nitermax aa b = do
   (q, h) <- arnoldi aa b nitermax
   (o, r) <- qr h
   return (q, o, extractDiagDense r)
@@ -372,7 +394,7 @@
 
 -- * Cholesky factorization
 
--- | Given a positive semidefinite matrix A, returns a lower-triangular matrix L such that L L^T = A . This is an implementation of the Cholesky–Banachiewicz algorithm, i.e. proceeding row by row from the upper-left corner.
+-- | Given a positive semidefinite matrix A, returns a lower-triangular matrix L such that L L^T = A . This is an implementation of the Cholesky–Banachiewicz algorithm, i.e. proceeding row by row from the upper-left corner and initializing the diagonal of the L matrix to ones.
 chol :: (Elt a, Epsilon a, MonadThrow m) =>
         SpMatrix a
      -> m (SpMatrix a)  -- ^ L
@@ -429,10 +451,8 @@
 
 
 -- * LU factorization
--- ** Doolittle algorithm
-{- Doolittle algorithm for factoring A' = P A, where P is a permutation matrix such that A' has a nonzero as its (0, 0) entry -}
 
--- | Given a matrix A, returns a pair of matrices (L, U) where L is lower triangular and U is upper triangular such that L U = A
+-- | Given a matrix A, returns a pair of matrices (L, U) where L is lower triangular and U is upper triangular such that L U = A . Implements the Doolittle algorithm, which expects all diagonal entries of A to be nonzero. Apply pivoting (row or column permutation) otherwise.
 lu :: (Scalar (SpVector t) ~ t, Elt t, VectorSpace (SpVector t), Epsilon t,
         MonadThrow m) =>
      SpMatrix t
@@ -581,7 +601,7 @@
 
 -- ** Incomplete LU
 
--- | Used for Incomplete LU : remove entries in `m` corresponding to zero entries in `m2` (this is called ILU(0) in the preconditioner literature)
+-- | Used for Incomplete LU : remove entries in the output matrix corresponding to zero entries in the input matrix (this is called ILU(0) in the preconditioner literature)
 ilu0 :: (Scalar (SpVector t) ~ t, Elt t, VectorSpace (SpVector t),
       Epsilon t, MonadThrow m) =>
      SpMatrix t
@@ -870,7 +890,7 @@
 
 
 -- * Moore-Penrose pseudoinverse
--- | Least-squares approximation of a rectangular system of equaitons. Uses <\\> for the linear solve
+-- | Least-squares approximation of a rectangular system of equaitons. Uses (<\>) for the linear solve
 pinv :: (MatrixType v ~ SpMatrix a, LinearSystem v, Epsilon a,
          MonadThrow m, MonadIO m) =>
      SpMatrix a -> v -> m v