anchor.diff:              for k=i to m-1 do
                   A.[l,k] <- A.[l,k] - 2.0 * v.[l-i] * v'A.[k-i]
           v                                                              // Return reflection vector.
  -            
  -    let QR (A:matrix) =
  +    
  +    type SparseEntry (index:int, value:float) =
  +        member m.index = index
  +        member val value = value with get, set
  +
  +    type SparseVector() = 
  +        member val entries = new ResizeArray<SparseEntry>() with get,set
  +        member m.Front =
  +            m.entries.[0]
  +        member m.Add (entry:SparseEntry) =
  +            m.entries.Add(entry)
  +        member m.IsEmpty =
  +            m.entries.Count <= 0
  +        member m.Swap (other:SparseVector) =
  +            let temp = new ResizeArray<SparseEntry> (other.entries)
  +            other.entries <- new ResizeArray<SparseEntry>(m.entries)
  +            m.entries <- new ResizeArray<SparseEntry>(temp)
  +        member m.PopFront =
  +            let value = m.entries.[0]
  +            m.entries.RemoveAt(0)
  +            value
  +
  +    /// QR factorization function for dense matrices
  +    let QRDense (A:matrix) =
           let (n,m) = matrixDims A
  -        let mutable Q = Matrix.identity n                                   // Keeps track of the orthogonal matrix.
  +        let mutable Q = Matrix.identity n // Keeps track of the orthogonal matrix.
           let R = Matrix.copy A
   
           // This method will update the orhogonal transformation fast when given a reflection vector.
  
positive.code:             UpdateQ Q v
          Q,R
  
      let Hessenberg (A:matrix) =
          // REVIEW we can do with less copying here.
          let (n,m) = matrixDims A
  
negative.code:                                         )
          leverage
          
  
      /// Calculates the pseudo inverse of the matrix
      let pseudoInvers (matrix:Matrix<float>) =
          let (m,n) = matrixDims matrix
  