anchor.diff:          (((*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
  
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:             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.
      
      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 R = Matrix.copy A
  
          // This method will update the orhogonal transformation fast when given a reflection vector.
  