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:         Q,Matrix.init n m (fun i j -> if i = 0 then A.[i,j] else R.[i-1,j])
          
      let leastSquares A (b: vector) =
          let (m,n) = matrixDims A
          // Is this an overdetermined or underdetermined system?
          if m > n then
              let qm,R = QR A
              let Qtb = qm.Transpose * b
              SolveTriangularLinearSystem R.[0..n-1,0..n-1] Qtb.[0..n-1] false
          else
              let qm,R = QR A
              let Qtb = qm.Transpose * b
              let s = SolveTriangularLinearSystem R.[0..m-1,0..m-1] Qtb false
              Vector.init n (fun i -> if i < m then s.[i] else 0.0)
  
      /// computes the hat matrix by the QR decomposition of the designmatrix used in ordinary least squares approaches
      let hatMatrix (designMatrix: Matrix<float>) = 
  
negative.code:         (((*P.Length,*)Permutation.ofArray P),L + Matrix.identity nA,U)
          //(P, L + (Matrix.identity nA), U)
      
      /// Solves a system of linear equations, AX = B, with A LU factorized.
      let SolveLinearSystem (A:matrix) (b:vector) =
          let (n,m) = matrixDims A
  