diff --git a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
index 979b362c..886c595c 100644
--- a/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
+++ b/src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.cs
@@ -229,7 +229,6 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             CommonParallel.For(0, y.Length, index => { result[index] = x[index] * y[index]; });
         }
 
-
         /// <summary>
         /// Does a point wise division of two arrays <c>z = x / y</c>. This can be used
         /// to divide elements of vectors or matrices.
@@ -289,8 +288,7 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 case Norm.FrobeniusNorm:
                     break;
             }
-
-            return ret;
+            throw new NotImplementedException();
         }
 
         /// <summary>
@@ -1206,36 +1204,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentNullException("a");
             }
 
-            for (var j = 0; j < order; j++)
+            var tmpColumn = new double[order];
+
+            // Main loop - along the diagonal
+            for (int ij = 0; ij < order; ij++)
             {
-                var d = 0.0;
-                int index;
-                for (var k = 0; k < j; k++)
+                // "Pivot" element
+                double tmpVal = a[(ij * order) + ij];
+
+                if (tmpVal > 0.0)
                 {
-                    var s = 0.0;
-                    int i;
-                    for (i = 0; i < k; i++)
+                    tmpVal = Math.Sqrt(tmpVal);
+                    a[(ij * order) + ij] = tmpVal;
+                    tmpColumn[ij] = tmpVal;
+
+                    // Calculate multipliers and copy to local column
+                    // Current column, below the diagonal
+                    for (int i = ij + 1; i < order; i++)
                     {
-                        s += a[(i * order) + k] * a[(i * order) + j];
+                        a[(ij * order) + i] /= tmpVal;
+                        tmpColumn[i] = a[(ij * order) + i];
                     }
 
-                    var tmp = k * order;
-                    index = tmp + j;
-                    a[index] = s = (a[index] - s) / a[tmp + k];
-                    d += s * s;
+                    // Remaining columns, below the diagonal
+                    DoCholeskyStep(a, order, ij + 1, order, tmpColumn, Control.NumberOfParallelWorkerThreads);
                 }
-
-                index = (j * order) + j;
-                d = a[index] - d;
-                if (d <= 0.0)
+                else
                 {
                     throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
                 }
 
-                a[index] = Math.Sqrt(d);
-                for (var k = j + 1; k < order; k++)
+                for (int i = ij + 1; i < order; i++)
                 {
-                    a[(k * order) + j] = 0.0;
+                    a[(i * order) + ij] = 0.0;
+                }
+            }
+        }
+
+        /// <summary>
+        /// Calculate Cholesky step
+        /// </summary>
+        /// <param name="data">Factor matrix</param>
+        /// <param name="rowDim">Number of rows</param>
+        /// <param name="firstCol">Column start</param>
+        /// <param name="colLimit">Total columns</param>
+        /// <param name="multipliers">Multipliears calculated previously</param>
+        /// <param name="availableCores">Number of available processors</param>
+        private static void DoCholeskyStep(double[] data, int rowDim, int firstCol, int colLimit, double[] multipliers, int availableCores)
+        {
+            var tmpColCount = colLimit - firstCol;
+
+            if ((availableCores > 1) && (tmpColCount > 200))
+            {
+                var tmpSplit = firstCol + (tmpColCount / 3);
+                var tmpCores = availableCores / 2;
+
+                CommonParallel.Invoke(
+                    () => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
+                    () => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
+            }
+            else
+            {
+                for (var j = firstCol; j < colLimit; j++)
+                {
+                    var tmpVal = multipliers[j];
+                    for (var i = j; i < rowDim; i++)
+                    {
+                        data[(j * rowDim) + i] -= multipliers[i] * tmpVal;
+                    }
                 }
             }
         }
@@ -1428,13 +1464,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             var minmn = Math.Min(rowsR, columnsR);
             for (var i = 0; i < minmn; i++)
             {
-                GenerateColumn(work, r, rowsR, i, rowsR - 1, i);
-                ComputeQR(work, i, r, rowsR, i, rowsR - 1, i + 1, columnsR - 1);
+                GenerateColumn(work, r, rowsR, i, i);
+                ComputeQR(work, i, r, i, rowsR, i + 1, columnsR, Control.NumberOfParallelWorkerThreads);
             }
 
             for (var i = minmn - 1; i >= 0; i--)
             {
-                ComputeQR(work, i, q, rowsR, i, rowsR - 1, i, rowsR - 1);
+                ComputeQR(work, i, q, i, rowsR, i, rowsR, Control.NumberOfParallelWorkerThreads);
             }
 
             work[0] = rowsR * rowsR;
@@ -1448,32 +1484,43 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="work">Work array</param>
         /// <param name="workIndex">Index of colunn in work array</param>
         /// <param name="a">Q or R matrices</param>
-        /// <param name="rowCount">The number of rows</param>
         /// <param name="rowStart">The first row in </param>
-        /// <param name="rowEnd">The last row</param>
+        /// <param name="rowCount">The last row</param>
         /// <param name="columnStart">The first column</param>
-        /// <param name="columnEnd">The last column</param>
-        private static void ComputeQR(double[] work, int workIndex, double[] a, int rowCount, int rowStart, int rowEnd, int columnStart, int columnEnd)
+        /// <param name="columnCount">The last column</param>
+        /// <param name="availableCores">Number of available CPUs</param>
+        private static void ComputeQR(double[] work, int workIndex, double[] a, int rowStart, int rowCount, int columnStart, int columnCount, int availableCores)
         {
-            if (rowStart > rowEnd || columnStart > columnEnd)
+            if (rowStart > rowCount || columnStart > columnCount)
             {
                 return;
             }
 
-            var vector = new double[columnEnd - columnStart + 1];
-            for (var i = rowStart; i <= rowEnd; i++)
+            var tmpColCount = columnCount - columnStart;
+
+            if ((availableCores > 1) && (tmpColCount > 200))
             {
-                for (var j = columnStart; j <= columnEnd; j++)
-                {
-                    vector[j - columnStart] += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
-                }
-            }
+                var tmpSplit = columnStart + (tmpColCount / 2);
+                var tmpCores = availableCores / 2;
 
-            for (var i = rowStart; i <= rowEnd; i++)
+                CommonParallel.Invoke(
+                    () => ComputeQR(work, workIndex, a, rowStart, rowCount, columnStart, tmpSplit, tmpCores),
+                    () => ComputeQR(work, workIndex, a, rowStart, rowCount, tmpSplit, columnCount, tmpCores));
+            }
+            else
             {
-                for (var j = columnStart; j <= columnEnd; j++)
+                for (var j = columnStart; j < columnCount; j++)
                 {
-                    a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart] * vector[j - columnStart];
+                    var scale = 0.0;
+                    for (var i = rowStart; i < rowCount; i++)     
+                    {
+                        scale += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
+                    }
+                
+                    for (var i = rowStart; i < rowCount; i++)    
+                    {
+                        a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart] * scale;
+                    }
                 }
             }
         }
@@ -1484,33 +1531,32 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="work">Work array</param>
         /// <param name="a">Initial matrix</param>
         /// <param name="rowCount">The number of rows in matrix</param>
-        /// <param name="rowStart">The firts row</param>
-        /// <param name="rowEnd">The last row</param>
+        /// <param name="row">The firts row</param>
         /// <param name="column">Column index</param>
-        private static void GenerateColumn(double[] work, double[] a, int rowCount, int rowStart, int rowEnd, int column)
+        private static void GenerateColumn(double[] work, double[] a, int rowCount, int row, int column)
         {
             var tmp = column * rowCount;
-            var index = tmp + rowStart;
+            var index = tmp + row;
 
             CommonParallel.For(
-                rowStart,
-                rowEnd + 1,
+                row,
+                rowCount,
                 i =>
                 {
                     var iIndex = tmp + i;
-                    work[iIndex - rowStart] = a[iIndex];
+                    work[iIndex - row] = a[iIndex];
                     a[iIndex] = 0.0;
                 });
 
             var norm = 0.0;
-            for (var i = 0; i < rowEnd - rowStart + 1; ++i)
+            for (var i = 0; i < rowCount - row; ++i)
             {
                 var iindex = tmp + i;
                 norm += work[iindex] * work[iindex];
             }
 
             norm = Math.Sqrt(norm);
-            if (rowStart == rowEnd || norm == 0)
+            if (row == rowCount - 1  || norm == 0)
             {
                 a[index] = -work[tmp];
                 work[tmp] = Math.Sqrt(2.0);
@@ -1524,11 +1570,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             }
 
             a[index] = -1.0 / scale;
-            CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] *= scale);
+            CommonParallel.For(0, rowCount - row, i => work[tmp + i] *= scale);
             work[tmp] += 1.0;
 
             var s = Math.Sqrt(1.0 / work[tmp]);
-            CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] *= s);
+            CommonParallel.For(0, rowCount - row, i => work[tmp + i] *= s);
         }
 
         #endregion
@@ -3977,43 +4023,81 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <remarks>This is equivalent to the POTRF LAPACK routine.</remarks>
         public void CholeskyFactor(float[] a, int order)
         {
-            var factor = new float[a.Length];
+            if (a == null)
+            {
+                throw new ArgumentNullException("a");
+            }
 
-            for (var j = 0; j < order; j++)
+            var tmpColumn = new float[order];
+
+            // Main loop - along the diagonal
+            for (var ij = 0; ij < order; ij++)
             {
-                float d = 0.0f;
-                int index;
-                for (var k = 0; k < j; k++)
+                // "Pivot" element
+                var tmpVal = a[(ij * order) + ij];
+
+                if (tmpVal > 0.0)
                 {
-                    float s = 0.0f;
-                    int i;
-                    for (i = 0; i < k; i++)
+                    tmpVal = (float)Math.Sqrt(tmpVal);
+                    a[(ij * order) + ij] = tmpVal;
+                    tmpColumn[ij] = tmpVal;
+
+                    // Calculate multipliers and copy to local column
+                    // Current column, below the diagonal
+                    for (var i = ij + 1; i < order; i++)
                     {
-                        s += factor[(i * order) + k] * factor[(i * order) + j];
+                        a[(ij * order) + i] /= tmpVal;
+                        tmpColumn[i] = a[(ij * order) + i];
                     }
 
-                    var tmp = k * order;
-                    index = tmp + j;
-                    s = (a[index] - s) / factor[tmp + k];
-                    factor[index] = s;
-                    d += s * s;
+                    // Remaining columns, below the diagonal
+                    DoCholeskyStep(a, order, ij + 1, order, tmpColumn, Control.NumberOfParallelWorkerThreads);
                 }
-
-                index = (j * order) + j;
-                d = a[index] - d;
-                if (d <= 0.0F)
+                else
                 {
                     throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
                 }
 
-                factor[index] = (float)Math.Sqrt(d);
-                for (var k = j + 1; k < order; k++)
+                for (int i = ij + 1; i < order; i++)
                 {
-                    factor[(k * order) + j] = 0.0F;
+                    a[(i * order) + ij] = 0.0f;
                 }
             }
+        }
+
+        /// <summary>
+        /// Calculate Cholesky step
+        /// </summary>
+        /// <param name="data">Factor matrix</param>
+        /// <param name="rowDim">Number of rows</param>
+        /// <param name="firstCol">Column start</param>
+        /// <param name="colLimit">Total columns</param>
+        /// <param name="multipliers">Multipliears calculated previously</param>
+        /// <param name="availableCores">Number of available processors</param>
+        private static void DoCholeskyStep(float[] data, int rowDim, int firstCol, int colLimit, float[] multipliers, int availableCores)
+        {
+            var tmpColCount = colLimit - firstCol;
+
+            if ((availableCores > 1) && (tmpColCount > 200))
+            {
+                var tmpSplit = firstCol + (tmpColCount / 3);
+                var tmpCores = availableCores / 2;
 
-            Buffer.BlockCopy(factor, 0, a, 0, factor.Length * Constants.SizeOfFloat);
+                CommonParallel.Invoke(
+                    () => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
+                    () => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
+            }
+            else
+            {
+                for (var j = firstCol; j < colLimit; j++)
+                {
+                    var tmpVal = multipliers[j];
+                    for (var i = j; i < rowDim; i++)
+                    {
+                        data[(j * rowDim) + i] -= multipliers[i] * tmpVal;
+                    }
+                }
+            }
         }
 
         /// <summary>
@@ -4204,13 +4288,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             var minmn = Math.Min(rowsR, columnsR);
             for (var i = 0; i < minmn; i++)
             {
-                GenerateColumn(work, r, rowsR, i, rowsR - 1, i);
-                ComputeQR(work, i, r, rowsR, i, rowsR - 1, i + 1, columnsR - 1);
+                GenerateColumn(work, r, rowsR, i, i);
+                ComputeQR(work, i, r, i, rowsR, i + 1, columnsR, Control.NumberOfParallelWorkerThreads);
             }
 
             for (var i = minmn - 1; i >= 0; i--)
             {
-                ComputeQR(work, i, q, rowsR, i, rowsR - 1, i, rowsR - 1);
+                ComputeQR(work, i, q, i, rowsR, i, rowsR, Control.NumberOfParallelWorkerThreads);
             }
 
             work[0] = rowsR * rowsR;
@@ -4224,32 +4308,43 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="work">Work array</param>
         /// <param name="workIndex">Index of colunn in work array</param>
         /// <param name="a">Q or R matrices</param>
-        /// <param name="rowCount">The number of rows</param>
         /// <param name="rowStart">The first row in </param>
-        /// <param name="rowEnd">The last row</param>
+        /// <param name="rowCount">The last row</param>
         /// <param name="columnStart">The first column</param>
-        /// <param name="columnEnd">The last column</param>
-        private static void ComputeQR(float[] work, int workIndex, float[] a, int rowCount, int rowStart, int rowEnd, int columnStart, int columnEnd)
+        /// <param name="columnCount">The last column</param>
+        /// <param name="availableCores">Number of available CPUs</param>
+        private static void ComputeQR(float[] work, int workIndex, float[] a, int rowStart, int rowCount, int columnStart, int columnCount, int availableCores)
         {
-            if (rowStart > rowEnd || columnStart > columnEnd)
+            if (rowStart > rowCount || columnStart > columnCount)
             {
                 return;
             }
 
-            var vector = new float[columnEnd - columnStart + 1];
-            for (var i = rowStart; i <= rowEnd; i++)
+            var tmpColCount = columnCount - columnStart;
+
+            if ((availableCores > 1) && (tmpColCount > 200))
             {
-                for (var j = columnStart; j <= columnEnd; j++)
-                {
-                    vector[j - columnStart] += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
-                }
-            }
+                var tmpSplit = columnStart + (tmpColCount / 2);
+                var tmpCores = availableCores / 2;
 
-            for (var i = rowStart; i <= rowEnd; i++)
+                CommonParallel.Invoke(
+                    () => ComputeQR(work, workIndex, a, rowStart, rowCount, columnStart, tmpSplit, tmpCores),
+                    () => ComputeQR(work, workIndex, a, rowStart, rowCount, tmpSplit, columnCount, tmpCores));
+            }
+            else
             {
-                for (var j = columnStart; j <= columnEnd; j++)
+                for (var j = columnStart; j < columnCount; j++)
                 {
-                    a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart] * vector[j - columnStart];
+                    var scale = 0.0f;
+                    for (var i = rowStart; i < rowCount; i++)
+                    {
+                        scale += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
+                    }
+
+                    for (var i = rowStart; i < rowCount; i++)
+                    {
+                        a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart] * scale;
+                    }
                 }
             }
         }
@@ -4260,33 +4355,32 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="work">Work array</param>
         /// <param name="a">Initial matrix</param>
         /// <param name="rowCount">The number of rows in matrix</param>
-        /// <param name="rowStart">The firts row</param>
-        /// <param name="rowEnd">The last row</param>
+        /// <param name="row">The firts row</param>
         /// <param name="column">Column index</param>
-        private static void GenerateColumn(float[] work, float[] a, int rowCount, int rowStart, int rowEnd, int column)
+        private static void GenerateColumn(float[] work, float[] a, int rowCount, int row, int column)
         {
             var tmp = column * rowCount;
-            var index = tmp + rowStart;
+            var index = tmp + row;
 
             CommonParallel.For(
-                rowStart,
-                rowEnd + 1,
+                row,
+                rowCount,
                 i =>
                 {
                     var iIndex = tmp + i;
-                    work[iIndex - rowStart] = a[iIndex];
+                    work[iIndex - row] = a[iIndex];
                     a[iIndex] = 0.0f;
                 });
 
             var norm = 0.0;
-            for (var i = 0; i < rowEnd - rowStart + 1; ++i)
+            for (var i = 0; i < rowCount - row; ++i)
             {
                 var iindex = tmp + i;
                 norm += work[iindex] * work[iindex];
             }
 
             norm = Math.Sqrt(norm);
-            if (rowStart == rowEnd || norm == 0)
+            if (row == rowCount - 1 || norm == 0)
             {
                 a[index] = -work[tmp];
                 work[tmp] = (float)Math.Sqrt(2.0);
@@ -4300,11 +4394,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             }
 
             a[index] = -1.0f / scale;
-            CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] *= scale);
+            CommonParallel.For(0, rowCount - row, i => work[tmp + i] *= scale);
             work[tmp] += 1.0f;
 
             var s = (float)Math.Sqrt(1.0 / work[tmp]);
-            CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] *= s);
+            CommonParallel.For(0, rowCount - row, i => work[tmp + i] *= s);
         }
 
         #endregion
@@ -6769,35 +6863,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentNullException("a");
             }
 
-            for (var i = 0; i < order; i++)
+            var tmpColumn = new Complex[order];
+
+            // Main loop - along the diagonal
+            for (var ij = 0; ij < order; ij++)
             {
-                var d = Complex.Zero;
-                int index;
-                for (var j = 0; j < i; j++)
+                // "Pivot" element
+                var tmpVal = a[(ij * order) + ij];
+
+                if (tmpVal.Real > 0.0)
                 {
-                    var s = Complex.Zero;
-                    for (var k = 0; k < j; k++)
+                    tmpVal = tmpVal.SquareRoot();
+                    a[(ij * order) + ij] = tmpVal;
+                    tmpColumn[ij] = tmpVal;
+
+                    // Calculate multipliers and copy to local column
+                    // Current column, below the diagonal
+                    for (var i = ij + 1; i < order; i++)
                     {
-                        s += a[(k * order) + i] * a[(k * order) + j].Conjugate();
+                        a[(ij * order) + i] /= tmpVal;
+                        tmpColumn[i] = a[(ij * order) + i];
                     }
 
-                    var tmp = j * order;
-                    index = tmp + i;
-                    a[index] = s = (a[index] - s) / a[tmp + j];
-                    d += s * s.Conjugate();
+                    // Remaining columns, below the diagonal
+                    DoCholeskyStep(a, order, ij + 1, order, tmpColumn, Control.NumberOfParallelWorkerThreads);
                 }
-
-                index = (i * order) + i;
-                d = a[index] - d;
-                if (d.Real <= 0.0)
+                else
                 {
                     throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
                 }
 
-                a[index] = d.SquareRoot();
-                for (var k = i + 1; k < order; k++)
+                for (var i = ij + 1; i < order; i++)
                 {
-                    a[(k * order) + i] = 0.0;
+                    a[(i * order) + ij] = 0.0;
+                }
+            }
+        }
+
+        /// <summary>
+        /// Calculate Cholesky step
+        /// </summary>
+        /// <param name="data">Factor matrix</param>
+        /// <param name="rowDim">Number of rows</param>
+        /// <param name="firstCol">Column start</param>
+        /// <param name="colLimit">Total columns</param>
+        /// <param name="multipliers">Multipliears calculated previously</param>
+        /// <param name="availableCores">Number of available processors</param>
+        private static void DoCholeskyStep(Complex[] data, int rowDim, int firstCol, int colLimit, Complex[] multipliers, int availableCores)
+        {
+            var tmpColCount = colLimit - firstCol;
+
+            if ((availableCores > 1) && (tmpColCount > 200))
+            {
+                var tmpSplit = firstCol + (tmpColCount / 3);
+                var tmpCores = availableCores / 2;
+
+                CommonParallel.Invoke(
+                    () => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
+                    () => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
+            }
+            else
+            {
+                for (var j = firstCol; j < colLimit; j++)
+                {
+                    var tmpVal = multipliers[j];
+                    for (var i = j; i < rowDim; i++)
+                    {
+                        data[(j * rowDim) + i] -= multipliers[i] * tmpVal.Conjugate();
+                    }
                 }
             }
         }
@@ -6990,13 +7123,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             var minmn = Math.Min(rowsR, columnsR);
             for (var i = 0; i < minmn; i++)
             {
-                GenerateColumn(work, r, rowsR, i, rowsR - 1, i);
-                ComputeQR(work, i, r, rowsR, i, rowsR - 1, i + 1, columnsR - 1);
+                GenerateColumn(work, r, rowsR, i, i);
+                ComputeQR(work, i, r, i, rowsR, i + 1, columnsR, Control.NumberOfParallelWorkerThreads);
             }
 
             for (var i = minmn - 1; i >= 0; i--)
             {
-                ComputeQR(work, i, q, rowsR, i, rowsR - 1, i, rowsR - 1);
+                ComputeQR(work, i, q, i, rowsR, i, rowsR, Control.NumberOfParallelWorkerThreads);
             }
 
             work[0] = rowsR * rowsR;
@@ -7010,32 +7143,43 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="work">Work array</param>
         /// <param name="workIndex">Index of colunn in work array</param>
         /// <param name="a">Q or R matrices</param>
-        /// <param name="rowCount">The number of rows</param>
         /// <param name="rowStart">The first row in </param>
-        /// <param name="rowEnd">The last row</param>
+        /// <param name="rowCount">The last row</param>
         /// <param name="columnStart">The first column</param>
-        /// <param name="columnEnd">The last column</param>
-        private static void ComputeQR(Complex[] work, int workIndex, Complex[] a, int rowCount, int rowStart, int rowEnd, int columnStart, int columnEnd)
+        /// <param name="columnCount">The last column</param>
+        /// <param name="availableCores">Number of available CPUs</param>
+        private static void ComputeQR(Complex[] work, int workIndex, Complex[] a, int rowStart, int rowCount, int columnStart, int columnCount, int availableCores)
         {
-            if (rowStart > rowEnd || columnStart > columnEnd)
+            if (rowStart > rowCount || columnStart > columnCount)
             {
                 return;
             }
 
-            var vector = new Complex[columnEnd - columnStart + 1];
-            for (var i = rowStart; i <= rowEnd; i++)
+            var tmpColCount = columnCount - columnStart;
+
+            if ((availableCores > 1) && (tmpColCount > 200))
             {
-                for (var j = columnStart; j <= columnEnd; j++)
-                {
-                    vector[j - columnStart] += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
-                }
-            }
+                var tmpSplit = columnStart + (tmpColCount / 2);
+                var tmpCores = availableCores / 2;
 
-            for (var i = rowStart; i <= rowEnd; i++)
+                CommonParallel.Invoke(
+                    () => ComputeQR(work, workIndex, a, rowStart, rowCount, columnStart, tmpSplit, tmpCores),
+                    () => ComputeQR(work, workIndex, a, rowStart, rowCount, tmpSplit, columnCount, tmpCores));
+            }
+            else
             {
-                for (var j = columnStart; j <= columnEnd; j++)
+                for (var j = columnStart; j < columnCount; j++)
                 {
-                    a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart].Conjugate() * vector[j - columnStart];
+                    var scale = Complex.Zero;
+                    for (var i = rowStart; i < rowCount; i++)
+                    {
+                        scale += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
+                    }
+
+                    for (var i = rowStart; i < rowCount; i++)
+                    {
+                        a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart].Conjugate() * scale;
+                    }
                 }
             }
         }
@@ -7046,33 +7190,32 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="work">Work array</param>
         /// <param name="a">Initial matrix</param>
         /// <param name="rowCount">The number of rows in matrix</param>
-        /// <param name="rowStart">The firts row</param>
-        /// <param name="rowEnd">The last row</param>
+        /// <param name="row">The firts row</param>
         /// <param name="column">Column index</param>
-        private static void GenerateColumn(Complex[] work, Complex[] a, int rowCount, int rowStart, int rowEnd, int column)
+        private static void GenerateColumn(Complex[] work, Complex[] a, int rowCount, int row, int column)
         {
             var tmp = column * rowCount;
-            var index = tmp + rowStart;
+            var index = tmp + row;
 
             CommonParallel.For(
-                rowStart,
-                rowEnd + 1,
+                row,
+                rowCount,
                 i =>
                 {
                     var iIndex = tmp + i;
-                    work[iIndex - rowStart] = a[iIndex];
+                    work[iIndex - row] = a[iIndex];
                     a[iIndex] = Complex.Zero;
                 });
 
             var norm = Complex.Zero;
-            for (var i = 0; i < rowEnd - rowStart + 1; ++i)
+            for (var i = 0; i < rowCount - row; ++i)
             {
                 var index1 = tmp + i;
                 norm += work[index1].Magnitude * work[index1].Magnitude;
             }
 
             norm = norm.SquareRoot();
-            if (rowStart == rowEnd || norm.Magnitude == 0)
+            if (row == rowCount - 1 || norm.Magnitude == 0)
             {
                 a[index] = -work[tmp];
                 work[tmp] = new Complex(2.0, 0).SquareRoot();
@@ -7085,11 +7228,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             }
 
             a[index] = -norm;
-            CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] /= norm);
+            CommonParallel.For(0, rowCount - row, i => work[tmp + i] /= norm);
             work[tmp] += 1.0;
 
             var s = (1.0 / work[tmp]).SquareRoot();
-            CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] = work[tmp + i].Conjugate() * s);
+            CommonParallel.For(0, rowCount - row, i => work[tmp + i] = work[tmp + i].Conjugate() * s);
         }
 
         #endregion
@@ -9505,35 +9648,74 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
                 throw new ArgumentNullException("a");
             }
 
-            for (var i = 0; i < order; i++)
+            var tmpColumn = new Complex32[order];
+
+            // Main loop - along the diagonal
+            for (var ij = 0; ij < order; ij++)
             {
-                var d = Complex32.Zero;
-                int index;
-                for (var j = 0; j < i; j++)
+                // "Pivot" element
+                var tmpVal = a[(ij * order) + ij];
+
+                if (tmpVal.Real > 0.0)
                 {
-                    var s = Complex32.Zero;
-                    for (var k = 0; k < j; k++)
+                    tmpVal = tmpVal.SquareRoot();
+                    a[(ij * order) + ij] = tmpVal;
+                    tmpColumn[ij] = tmpVal;
+
+                    // Calculate multipliers and copy to local column
+                    // Current column, below the diagonal
+                    for (var i = ij + 1; i < order; i++)
                     {
-                        s += a[(k * order) + i] * a[(k * order) + j].Conjugate();
+                        a[(ij * order) + i] /= tmpVal;
+                        tmpColumn[i] = a[(ij * order) + i];
                     }
 
-                    var tmp = j * order;
-                    index = tmp + i;
-                    a[index] = s = (a[index] - s) / a[tmp + j];
-                    d += s * s.Conjugate();
+                    // Remaining columns, below the diagonal
+                    DoCholeskyStep(a, order, ij + 1, order, tmpColumn, Control.NumberOfParallelWorkerThreads);
                 }
-
-                index = (i * order) + i;
-                d = a[index] - d;
-                if (d.Real <= 0.0f)
+                else
                 {
                     throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
                 }
 
-                a[index] = d.SquareRoot();
-                for (var k = i + 1; k < order; k++)
+                for (var i = ij + 1; i < order; i++)
                 {
-                    a[(k * order) + i] = 0.0f;
+                    a[(i * order) + ij] = 0.0f;
+                }
+            }
+        }
+
+        /// <summary>
+        /// Calculate Cholesky step
+        /// </summary>
+        /// <param name="data">Factor matrix</param>
+        /// <param name="rowDim">Number of rows</param>
+        /// <param name="firstCol">Column start</param>
+        /// <param name="colLimit">Total columns</param>
+        /// <param name="multipliers">Multipliears calculated previously</param>
+        /// <param name="availableCores">Number of available processors</param>
+        private static void DoCholeskyStep(Complex32[] data, int rowDim, int firstCol, int colLimit, Complex32[] multipliers, int availableCores)
+        {
+            var tmpColCount = colLimit - firstCol;
+
+            if ((availableCores > 1) && (tmpColCount > 200))
+            {
+                var tmpSplit = firstCol + (tmpColCount / 3);
+                var tmpCores = availableCores / 2;
+
+                CommonParallel.Invoke(
+                    () => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
+                    () => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
+            }
+            else
+            {
+                for (var j = firstCol; j < colLimit; j++)
+                {
+                    var tmpVal = multipliers[j];
+                    for (var i = j; i < rowDim; i++)
+                    {
+                        data[(j * rowDim) + i] -= multipliers[i] * tmpVal.Conjugate();
+                    }
                 }
             }
         }
@@ -9726,13 +9908,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             var minmn = Math.Min(rowsR, columnsR);
             for (var i = 0; i < minmn; i++)
             {
-                GenerateColumn(work, r, rowsR, i, rowsR - 1, i);
-                ComputeQR(work, i, r, rowsR, i, rowsR - 1, i + 1, columnsR - 1);
+                GenerateColumn(work, r, rowsR, i, i);
+                ComputeQR(work, i, r, i, rowsR, i + 1, columnsR, Control.NumberOfParallelWorkerThreads);
             }
 
             for (var i = minmn - 1; i >= 0; i--)
             {
-                ComputeQR(work, i, q, rowsR, i, rowsR - 1, i, rowsR - 1);
+                ComputeQR(work, i, q, i, rowsR, i, rowsR, Control.NumberOfParallelWorkerThreads);
             }
 
             work[0] = rowsR * rowsR;
@@ -9746,32 +9928,43 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="work">Work array</param>
         /// <param name="workIndex">Index of colunn in work array</param>
         /// <param name="a">Q or R matrices</param>
-        /// <param name="rowCount">The number of rows</param>
         /// <param name="rowStart">The first row in </param>
-        /// <param name="rowEnd">The last row</param>
+        /// <param name="rowCount">The last row</param>
         /// <param name="columnStart">The first column</param>
-        /// <param name="columnEnd">The last column</param>
-        private static void ComputeQR(Complex32[] work, int workIndex, Complex32[] a, int rowCount, int rowStart, int rowEnd, int columnStart, int columnEnd)
+        /// <param name="columnCount">The last column</param>
+        /// <param name="availableCores">Number of available CPUs</param>
+        private static void ComputeQR(Complex32[] work, int workIndex, Complex32[] a, int rowStart, int rowCount, int columnStart, int columnCount, int availableCores)
         {
-            if (rowStart > rowEnd || columnStart > columnEnd)
+            if (rowStart > rowCount || columnStart > columnCount)
             {
                 return;
             }
 
-            var vector = new Complex32[columnEnd - columnStart + 1];
-            for (var i = rowStart; i <= rowEnd; i++)
+            var tmpColCount = columnCount - columnStart;
+
+            if ((availableCores > 1) && (tmpColCount > 200))
             {
-                for (var j = columnStart; j <= columnEnd; j++)
-                {
-                    vector[j - columnStart] += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
-                }
-            }
+                var tmpSplit = columnStart + (tmpColCount / 2);
+                var tmpCores = availableCores / 2;
 
-            for (var i = rowStart; i <= rowEnd; i++)
+                CommonParallel.Invoke(
+                    () => ComputeQR(work, workIndex, a, rowStart, rowCount, columnStart, tmpSplit, tmpCores),
+                    () => ComputeQR(work, workIndex, a, rowStart, rowCount, tmpSplit, columnCount, tmpCores));
+            }
+            else
             {
-                for (var j = columnStart; j <= columnEnd; j++)
+                for (var j = columnStart; j < columnCount; j++)
                 {
-                    a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart].Conjugate() * vector[j - columnStart];
+                    var scale = Complex32.Zero;
+                    for (var i = rowStart; i < rowCount; i++)
+                    {
+                        scale += work[(workIndex * rowCount) + i - rowStart] * a[(j * rowCount) + i];
+                    }
+
+                    for (var i = rowStart; i < rowCount; i++)
+                    {
+                        a[(j * rowCount) + i] -= work[(workIndex * rowCount) + i - rowStart].Conjugate() * scale;
+                    }
                 }
             }
         }
@@ -9782,33 +9975,32 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="work">Work array</param>
         /// <param name="a">Initial matrix</param>
         /// <param name="rowCount">The number of rows in matrix</param>
-        /// <param name="rowStart">The firts row</param>
-        /// <param name="rowEnd">The last row</param>
+        /// <param name="row">The firts row</param>
         /// <param name="column">Column index</param>
-        private static void GenerateColumn(Complex32[] work, Complex32[] a, int rowCount, int rowStart, int rowEnd, int column)
+        private static void GenerateColumn(Complex32[] work, Complex32[] a, int rowCount, int row, int column)
         {
             var tmp = column * rowCount;
-            var index = tmp + rowStart;
+            var index = tmp + row;
 
             CommonParallel.For(
-                rowStart,
-                rowEnd + 1,
+                row,
+                rowCount,
                 i =>
                 {
                     var iIndex = tmp + i;
-                    work[iIndex - rowStart] = a[iIndex];
+                    work[iIndex - row] = a[iIndex];
                     a[iIndex] = Complex32.Zero;
                 });
 
             var norm = Complex32.Zero;
-            for (var i = 0; i < rowEnd - rowStart + 1; ++i)
+            for (var i = 0; i < rowCount - row; ++i)
             {
                 var index1 = tmp + i;
                 norm += work[index1].Magnitude * work[index1].Magnitude;
             }
 
             norm = norm.SquareRoot();
-            if (rowStart == rowEnd || norm.Magnitude == 0)
+            if (row == rowCount - 1 || norm.Magnitude == 0)
             {
                 a[index] = -work[tmp];
                 work[tmp] = new Complex32(2.0f, 0).SquareRoot();
@@ -9821,11 +10013,11 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
             }
 
             a[index] = -norm;
-            CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] /= norm);
+            CommonParallel.For(0, rowCount - row, i => work[tmp + i] /= norm);
             work[tmp] += 1.0f;
 
             var s = (1.0f / work[tmp]).SquareRoot();
-            CommonParallel.For(0, rowEnd - rowStart + 1, i => work[tmp + i] = work[tmp + i].Conjugate() * s);
+            CommonParallel.For(0, rowCount - row, i => work[tmp + i] = work[tmp + i].Conjugate() * s);
         }
 
         #endregion
diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs
index 56ea890f..ef82e2f0 100644
--- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs
+++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserCholesky.cs
@@ -33,8 +33,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
     using System;
     using System.Numerics;
     using Generic;
-    using Generic.Factorization;
     using Properties;
+    using Threading;
 
     /// <summary>
     /// <para>A class which encapsulates the functionality of a Cholesky factorization for user matrices.</para>
@@ -69,32 +69,74 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
 
             // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
             CholeskyFactor = matrix.Clone();
-            for (var i = 0; i < CholeskyFactor.RowCount; i++)
+            var tmpColumn = new Complex[CholeskyFactor.RowCount];
+
+            // Main loop - along the diagonal
+            for (var ij = 0; ij < CholeskyFactor.RowCount; ij++)
             {
-                var d = Complex.Zero;
-                for (var j = 0; j < i; j++)
+                // "Pivot" element
+                var tmpVal = CholeskyFactor.At(ij, ij);
+
+                if (tmpVal.Real > 0.0)
                 {
-                    var s = Complex.Zero;
-                    for (var k = 0; k < j; k++)
+                    tmpVal = tmpVal.SquareRoot();
+                    CholeskyFactor.At(ij, ij, tmpVal);
+                    tmpColumn[ij] = tmpVal;
+
+                    // Calculate multipliers and copy to local column
+                    // Current column, below the diagonal
+                    for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
                     {
-                        s += CholeskyFactor.At(i, k) * CholeskyFactor.At(j, k).Conjugate();
+                        CholeskyFactor.At(i, ij, CholeskyFactor.At(i, ij) / tmpVal);
+                        tmpColumn[i] = CholeskyFactor.At(i, ij);
                     }
 
-                    s = (matrix.At(i, j) - s) / CholeskyFactor.At(j, j);
-                    CholeskyFactor.At(i, j, s);
-                    d += s * s.Conjugate();
+                    // Remaining columns, below the diagonal
+                    DoCholeskyStep(CholeskyFactor, CholeskyFactor.RowCount, ij + 1, CholeskyFactor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads);
                 }
-
-                d = matrix.At(i, i) - d;
-                if (d.Real <= 0)
+                else
                 {
                     throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
                 }
 
-                CholeskyFactor.At(i, i, d.SquareRoot());
-                for (var k = i + 1; k < CholeskyFactor.RowCount; k++)
+                for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
                 {
-                    CholeskyFactor.At(i, k, 0.0);
+                    CholeskyFactor.At(ij, i, Complex.Zero);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Calculate Cholesky step
+        /// </summary>
+        /// <param name="data">Factor matrix</param>
+        /// <param name="rowDim">Number of rows</param>
+        /// <param name="firstCol">Column start</param>
+        /// <param name="colLimit">Total columns</param>
+        /// <param name="multipliers">Multipliers calculated previously</param>
+        /// <param name="availableCores">Number of available processors</param>
+        private static void DoCholeskyStep(Matrix<Complex> data, int rowDim, int firstCol, int colLimit, Complex[] multipliers, int availableCores)
+        {
+            var tmpColCount = colLimit - firstCol;
+
+            if ((availableCores > 1) && (tmpColCount > 200))
+            {
+                var tmpSplit = firstCol + (tmpColCount / 3);
+                var tmpCores = availableCores / 2;
+
+                CommonParallel.Invoke(
+                    () => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
+                    () => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
+            }
+            else
+            {
+                for (var j = firstCol; j < colLimit; j++)
+                {
+                    var tmpVal = multipliers[j];
+                    for (var i = j; i < rowDim; i++)
+                    {
+                        data.At(i, j, data.At(i, j) - (multipliers[i] * tmpVal.Conjugate()));
+                    }
                 }
             }
         }
diff --git a/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs
index e4c9a382..68db066d 100644
--- a/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs
+++ b/src/Numerics/LinearAlgebra/Complex/Factorization/UserQR.cs
@@ -35,6 +35,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
     using System.Numerics;
     using Generic;
     using Properties;
+    using Threading;
 
     /// <summary>
     /// <para>A class which encapsulates the functionality of the QR decomposition.</para>
@@ -77,13 +78,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
             var u = new Complex[minmn][];
             for (var i = 0; i < minmn; i++)
             {
-                u[i] = GenerateColumn(MatrixR, i, matrix.RowCount - 1, i);
-                ComputeQR(u[i], MatrixR, i, matrix.RowCount - 1, i + 1, matrix.ColumnCount - 1);
+                u[i] = GenerateColumn(MatrixR, i, i);
+                ComputeQR(u[i], MatrixR, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
             }
 
             for (var i = minmn - 1; i >= 0; i--)
             {
-                ComputeQR(u[i], MatrixQ, i, matrix.RowCount - 1, i, matrix.RowCount - 1);
+                ComputeQR(u[i], MatrixQ, i, matrix.RowCount, i, matrix.RowCount, Control.NumberOfParallelWorkerThreads);
             }
         }
 
@@ -91,27 +92,26 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         /// Generate column from initial matrix to work array
         /// </summary>
         /// <param name="a">Initial matrix</param>
-        /// <param name="rowStart">The first row</param>
-        /// <param name="rowEnd">The last row</param>
+        /// <param name="row">The first row</param>
         /// <param name="column">Column index</param>
         /// <returns>Generated vector</returns>
-        private static Complex[] GenerateColumn(Matrix<Complex> a, int rowStart, int rowEnd, int column)
+        private static Complex[] GenerateColumn(Matrix<Complex> a, int row, int column)
         {
-            var ru = rowEnd - rowStart + 1;
+            var ru = a.RowCount - row;
             var u = new Complex[ru];
 
-            for (var i = rowStart; i <= rowEnd; i++)
+            for (var i = row; i < a.RowCount; i++)
             {
-                u[i - rowStart] = a.At(i, column);
+                u[i - row] = a.At(i, column);
                 a.At(i, column, 0.0);
             }
 
             var norm = u.Aggregate(Complex.Zero, (current, t) => current + (t.Magnitude * t.Magnitude));
             norm = norm.SquareRoot();
 
-            if (rowStart == rowEnd || norm.Magnitude == 0)
+            if (row == a.RowCount - 1 || norm.Magnitude == 0)
             {
-                a.At(rowStart, column, -u[0]);
+                a.At(row, column, -u[0]);
                 u[0] = Math.Sqrt(2.0);
                 return u;
             }
@@ -121,7 +121,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
                 norm = norm.Magnitude * (u[0] / u[0].Magnitude);
             }
 
-            a.At(rowStart, column, -norm);
+            a.At(row, column, -norm);
 
             for (var i = 0; i < ru; i++)
             {
@@ -140,40 +140,47 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
         }
 
         /// <summary>
-        /// Compute matrix Q or R
+        /// Perform calculation of Q or R
         /// </summary>
-        /// <param name="u">H vectors values</param>
-        /// <param name="a">Matrix to calculate</param>
-        /// <param name="rowStart">Starting row index</param>
-        /// <param name="rowEnd">Ending row index</param>
-        /// <param name="columnStart">Starting column index</param>
-        /// <param name="columnEnd">Ending column index</param>
-        private static void ComputeQR(Complex[] u, Matrix<Complex> a, int rowStart, int rowEnd, int columnStart, int columnEnd)
+        /// <param name="u">Work array</param>
+        /// <param name="a">Q or R matrices</param>
+        /// <param name="rowStart">The first row</param>
+        /// <param name="rowDim">The last row</param>
+        /// <param name="columnStart">The first column</param>
+        /// <param name="columnDim">The last column</param>
+        /// <param name="availableCores">Number of available CPUs</param>
+        private static void ComputeQR(Complex[] u, Matrix<Complex> a, int rowStart, int rowDim, int columnStart, int columnDim, int availableCores)
         {
-            if (rowEnd < rowStart || columnEnd < columnStart)
+            if (rowDim < rowStart || columnDim < columnStart)
             {
                 return;
             }
 
-            var v = new Complex[columnEnd - columnStart + 1];
-            for (var j = columnStart; j <= columnEnd; j++)
-            {
-                v[j - columnStart] = 0.0;
-            }
+            var tmpColCount = columnDim - columnStart;
 
-            for (var i = rowStart; i <= rowEnd; i++)
+            if ((availableCores > 1) && (tmpColCount > 200))
             {
-                for (var j = columnStart; j <= columnEnd; j++)
-                {
-                    v[j - columnStart] = v[j - columnStart] + (u[i - rowStart] * a.At(i, j));
-                }
-            }
+                var tmpSplit = columnStart + (tmpColCount / 2);
+                var tmpCores = availableCores / 2;
 
-            for (var i = rowStart; i <= rowEnd; i++)
+                CommonParallel.Invoke(
+                    () => ComputeQR(u, a, rowStart, rowDim, columnStart, tmpSplit, tmpCores),
+                    () => ComputeQR(u, a, rowStart, rowDim, tmpSplit, columnDim, tmpCores));
+            }
+            else
             {
-                for (var j = columnStart; j <= columnEnd; j++)
+                for (var j = columnStart; j < columnDim; j++)
                 {
-                    a.At(i, j, a.At(i, j) - (u[i - rowStart].Conjugate() * v[j - columnStart]));
+                    var scale = Complex.Zero;
+                    for (var i = rowStart; i < rowDim; i++)
+                    {
+                        scale += u[i - rowStart] * a.At(i, j);
+                    }
+
+                    for (var i = rowStart; i < rowDim; i++)
+                    {
+                        a.At(i, j, a.At(i, j) - (u[i - rowStart].Conjugate() * scale));
+                    }
                 }
             }
         }
diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserCholesky.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserCholesky.cs
index f7414b41..c3866922 100644
--- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserCholesky.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserCholesky.cs
@@ -34,6 +34,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
     using Generic;
     using Numerics;
     using Properties;
+    using Threading;
 
     /// <summary>
     /// <para>A class which encapsulates the functionality of a Cholesky factorization for user matrices.</para>
@@ -68,32 +69,74 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
 
             // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
             CholeskyFactor = matrix.Clone();
-            for (var i = 0; i < CholeskyFactor.RowCount; i++)
+            var tmpColumn = new Complex32[CholeskyFactor.RowCount];
+
+            // Main loop - along the diagonal
+            for (var ij = 0; ij < CholeskyFactor.RowCount; ij++)
             {
-                var d = Complex32.Zero;
-                for (var j = 0; j < i; j++)
+                // "Pivot" element
+                var tmpVal = CholeskyFactor.At(ij, ij);
+
+                if (tmpVal.Real > 0.0)
                 {
-                    var s = Complex32.Zero;
-                    for (var k = 0; k < j; k++)
+                    tmpVal = tmpVal.SquareRoot();
+                    CholeskyFactor.At(ij, ij, tmpVal);
+                    tmpColumn[ij] = tmpVal;
+
+                    // Calculate multipliers and copy to local column
+                    // Current column, below the diagonal
+                    for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
                     {
-                        s += CholeskyFactor.At(i, k) * CholeskyFactor.At(j, k).Conjugate();
+                        CholeskyFactor.At(i, ij, CholeskyFactor.At(i, ij) / tmpVal);
+                        tmpColumn[i] = CholeskyFactor.At(i, ij);
                     }
 
-                    s = (matrix.At(i, j) - s) / CholeskyFactor.At(j, j);
-                    CholeskyFactor.At(i, j, s);
-                    d += s * s.Conjugate();
+                    // Remaining columns, below the diagonal
+                    DoCholeskyStep(CholeskyFactor, CholeskyFactor.RowCount, ij + 1, CholeskyFactor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads);
                 }
-
-                d = matrix.At(i, i) - d;
-                if (d.Real <= 0)
+                else
                 {
                     throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
                 }
 
-                CholeskyFactor.At(i, i, d.SquareRoot());
-                for (var k = i + 1; k < CholeskyFactor.RowCount; k++)
+                for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
                 {
-                    CholeskyFactor.At(i, k, 0.0f);
+                    CholeskyFactor.At(ij, i, Complex32.Zero);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Calculate Cholesky step
+        /// </summary>
+        /// <param name="data">Factor matrix</param>
+        /// <param name="rowDim">Number of rows</param>
+        /// <param name="firstCol">Column start</param>
+        /// <param name="colLimit">Total columns</param>
+        /// <param name="multipliers">Multipliers calculated previously</param>
+        /// <param name="availableCores">Number of available processors</param>
+        private static void DoCholeskyStep(Matrix<Complex32> data, int rowDim, int firstCol, int colLimit, Complex32[] multipliers, int availableCores)
+        {
+            var tmpColCount = colLimit - firstCol;
+
+            if ((availableCores > 1) && (tmpColCount > 200))
+            {
+                var tmpSplit = firstCol + (tmpColCount / 3);
+                var tmpCores = availableCores / 2;
+
+                CommonParallel.Invoke(
+                    () => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
+                    () => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
+            }
+            else
+            {
+                for (var j = firstCol; j < colLimit; j++)
+                {
+                    var tmpVal = multipliers[j];
+                    for (var i = j; i < rowDim; i++)
+                    {
+                        data.At(i, j, data.At(i, j) - (multipliers[i] * tmpVal.Conjugate()));
+                    }
                 }
             }
         }
diff --git a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserQR.cs
index 9d8b21ea..333c64b8 100644
--- a/src/Numerics/LinearAlgebra/Complex32/Factorization/UserQR.cs
+++ b/src/Numerics/LinearAlgebra/Complex32/Factorization/UserQR.cs
@@ -35,6 +35,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
     using Generic;
     using Numerics;
     using Properties;
+    using Threading;
 
     /// <summary>
     /// <para>A class which encapsulates the functionality of the QR decomposition.</para>
@@ -77,13 +78,13 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
             var u = new Complex32[minmn][];
             for (var i = 0; i < minmn; i++)
             {
-                u[i] = GenerateColumn(MatrixR, i, matrix.RowCount - 1, i);
-                ComputeQR(u[i], MatrixR, i, matrix.RowCount - 1, i + 1, matrix.ColumnCount - 1);
+                u[i] = GenerateColumn(MatrixR, i, i);
+                ComputeQR(u[i], MatrixR, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
             }
 
             for (var i = minmn - 1; i >= 0; i--)
             {
-                ComputeQR(u[i], MatrixQ, i, matrix.RowCount - 1, i, matrix.RowCount - 1);
+                ComputeQR(u[i], MatrixQ, i, matrix.RowCount, i, matrix.RowCount, Control.NumberOfParallelWorkerThreads);
             }
         }
 
@@ -91,27 +92,26 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         /// Generate column from initial matrix to work array
         /// </summary>
         /// <param name="a">Initial matrix</param>
-        /// <param name="rowStart">The first row</param>
-        /// <param name="rowEnd">The last row</param>
+        /// <param name="row">The first row</param>
         /// <param name="column">Column index</param>
         /// <returns>Generated vector</returns>
-        private static Complex32[] GenerateColumn(Matrix<Complex32> a, int rowStart, int rowEnd, int column)
+        private static Complex32[] GenerateColumn(Matrix<Complex32> a, int row, int column)
         {
-            var ru = rowEnd - rowStart + 1;
+            var ru = a.RowCount - row;
             var u = new Complex32[ru];
 
-            for (var i = rowStart; i <= rowEnd; i++)
+            for (var i = row; i < a.RowCount; i++)
             {
-                u[i - rowStart] = a.At(i, column);
+                u[i - row] = a.At(i, column);
                 a.At(i, column, 0.0f);
             }
 
             var norm = u.Aggregate(Complex32.Zero, (current, t) => current + (t.Magnitude * t.Magnitude));
             norm = norm.SquareRoot();
 
-            if (rowStart == rowEnd || norm.Magnitude == 0)
+            if (row == a.RowCount - 1 || norm.Magnitude == 0)
             {
-                a.At(rowStart, column, -u[0]);
+                a.At(row, column, -u[0]);
                 u[0] = (float)Math.Sqrt(2.0);
                 return u;
             }
@@ -121,7 +121,7 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
                 norm = norm.Magnitude * (u[0] / u[0].Magnitude);
             }
 
-            a.At(rowStart, column, -norm);
+            a.At(row, column, -norm);
 
             for (var i = 0; i < ru; i++)
             {
@@ -140,40 +140,47 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
         }
 
         /// <summary>
-        /// Compute matrix Q or R
+        /// Perform calculation of Q or R
         /// </summary>
-        /// <param name="u">H vectors values</param>
-        /// <param name="a">Matrix to calculate</param>
-        /// <param name="rowStart">Starting row index</param>
-        /// <param name="rowEnd">Ending row index</param>
-        /// <param name="columnStart">Starting column index</param>
-        /// <param name="columnEnd">Ending column index</param>
-        private static void ComputeQR(Complex32[] u, Matrix<Complex32> a, int rowStart, int rowEnd, int columnStart, int columnEnd)
+        /// <param name="u">Work array</param>
+        /// <param name="a">Q or R matrices</param>
+        /// <param name="rowStart">The first row</param>
+        /// <param name="rowDim">The last row</param>
+        /// <param name="columnStart">The first column</param>
+        /// <param name="columnDim">The last column</param>
+        /// <param name="availableCores">Number of available CPUs</param>
+        private static void ComputeQR(Complex32[] u, Matrix<Complex32> a, int rowStart, int rowDim, int columnStart, int columnDim, int availableCores)
         {
-            if (rowEnd < rowStart || columnEnd < columnStart)
+            if (rowDim < rowStart || columnDim < columnStart)
             {
                 return;
             }
 
-            var v = new Complex32[columnEnd - columnStart + 1];
-            for (var j = columnStart; j <= columnEnd; j++)
-            {
-                v[j - columnStart] = 0.0f;
-            }
+            var tmpColCount = columnDim - columnStart;
 
-            for (var i = rowStart; i <= rowEnd; i++)
+            if ((availableCores > 1) && (tmpColCount > 200))
             {
-                for (var j = columnStart; j <= columnEnd; j++)
-                {
-                    v[j - columnStart] = v[j - columnStart] + (u[i - rowStart] * a.At(i, j));
-                }
-            }
+                var tmpSplit = columnStart + (tmpColCount / 2);
+                var tmpCores = availableCores / 2;
 
-            for (var i = rowStart; i <= rowEnd; i++)
+                CommonParallel.Invoke(
+                    () => ComputeQR(u, a, rowStart, rowDim, columnStart, tmpSplit, tmpCores),
+                    () => ComputeQR(u, a, rowStart, rowDim, tmpSplit, columnDim, tmpCores));
+            }
+            else
             {
-                for (var j = columnStart; j <= columnEnd; j++)
+                for (var j = columnStart; j < columnDim; j++)
                 {
-                    a.At(i, j, a.At(i, j) - (u[i - rowStart].Conjugate() * v[j - columnStart]));
+                    var scale = Complex32.Zero;
+                    for (var i = rowStart; i < rowDim; i++)
+                    {
+                        scale += u[i - rowStart] * a.At(i, j);
+                    }
+
+                    for (var i = rowStart; i < rowDim; i++)
+                    {
+                        a.At(i, j, a.At(i, j) - (u[i - rowStart].Conjugate() * scale));
+                    }
                 }
             }
         }
diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs
index 28c719f1..a0ceebba 100644
--- a/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs
+++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserCholesky.cs
@@ -33,6 +33,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
     using System;
     using Generic;
     using Properties;
+    using Threading;
 
     /// <summary>
     /// <para>A class which encapsulates the functionality of a Cholesky factorization for user matrices.</para>
@@ -67,32 +68,74 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
 
             // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
             CholeskyFactor = matrix.Clone();
-            for (var j = 0; j < CholeskyFactor.RowCount; j++)
+            var tmpColumn = new double[CholeskyFactor.RowCount];
+
+            // Main loop - along the diagonal
+            for (var ij = 0; ij < CholeskyFactor.RowCount; ij++)
             {
-                var d = 0.0;
-                for (var k = 0; k < j; k++)
+                // "Pivot" element
+                var tmpVal = CholeskyFactor.At(ij, ij);
+
+                if (tmpVal > 0.0)
                 {
-                    var s = 0.0;
-                    for (var i = 0; i < k; i++)
+                    tmpVal = Math.Sqrt(tmpVal);
+                    CholeskyFactor.At(ij, ij, tmpVal);
+                    tmpColumn[ij] = tmpVal;
+
+                    // Calculate multipliers and copy to local column
+                    // Current column, below the diagonal
+                    for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
                     {
-                        s += CholeskyFactor.At(k, i) * CholeskyFactor.At(j, i);
+                        CholeskyFactor.At(i, ij, CholeskyFactor.At(i, ij) / tmpVal);
+                        tmpColumn[i] = CholeskyFactor.At(i, ij);
                     }
 
-                    s = (matrix.At(j, k) - s) / CholeskyFactor.At(k, k);
-                    CholeskyFactor.At(j, k, s);
-                    d += s * s;
+                    // Remaining columns, below the diagonal
+                    DoCholeskyStep(CholeskyFactor, CholeskyFactor.RowCount, ij + 1, CholeskyFactor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads);
                 }
-
-                d = matrix.At(j, j) - d;
-                if (d <= 0.0)
+                else
                 {
                     throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
                 }
 
-                CholeskyFactor.At(j, j, Math.Sqrt(d));
-                for (var k = j + 1; k < CholeskyFactor.RowCount; k++)
+                for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
                 {
-                    CholeskyFactor.At(j, k, 0.0);
+                    CholeskyFactor.At(ij, i, 0.0);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Calculate Cholesky step
+        /// </summary>
+        /// <param name="data">Factor matrix</param>
+        /// <param name="rowDim">Number of rows</param>
+        /// <param name="firstCol">Column start</param>
+        /// <param name="colLimit">Total columns</param>
+        /// <param name="multipliers">Multipliers calculated previously</param>
+        /// <param name="availableCores">Number of available processors</param>
+        private static void DoCholeskyStep(Matrix<double> data, int rowDim, int firstCol, int colLimit, double[] multipliers, int availableCores)
+        {
+            var tmpColCount = colLimit - firstCol;
+
+            if ((availableCores > 1) && (tmpColCount > 200))
+            {
+                var tmpSplit = firstCol + (tmpColCount / 3);
+                var tmpCores = availableCores / 2;
+
+                CommonParallel.Invoke(
+                    () => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
+                    () => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
+            }
+            else
+            {
+                for (var j = firstCol; j < colLimit; j++)
+                {
+                    var tmpVal = multipliers[j];
+                    for (var i = j; i < rowDim; i++)
+                    {
+                        data.At(i, j, data.At(i, j) - (multipliers[i] * tmpVal));
+                    }
                 }
             }
         }
diff --git a/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs
index c9d6695c..44d3c8d1 100644
--- a/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs
+++ b/src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs
@@ -34,6 +34,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
     using System.Linq;
     using Generic;
     using Properties;
+    using Threading;
 
     /// <summary>
     /// <para>A class which encapsulates the functionality of the QR decomposition.</para>
@@ -76,41 +77,40 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
             var u = new double[minmn][];
             for (var i = 0; i < minmn; i++)
             {
-                u[i] = GenerateColumn(MatrixR, i, matrix.RowCount - 1, i);
-                ComputeQR(u[i], MatrixR, i, matrix.RowCount - 1, i + 1, matrix.ColumnCount - 1);
+                u[i] = GenerateColumn(MatrixR, i, i);
+                ComputeQR(u[i], MatrixR, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
             }
 
             for (var i = minmn - 1; i >= 0; i--)
             {
-                ComputeQR(u[i], MatrixQ, i, matrix.RowCount - 1, i, matrix.RowCount - 1);
+                ComputeQR(u[i], MatrixQ, i, matrix.RowCount, i, matrix.RowCount, Control.NumberOfParallelWorkerThreads);
             }
         }
-
+        
         /// <summary>
         /// Generate column from initial matrix to work array
         /// </summary>
         /// <param name="a">Initial matrix</param>
-        /// <param name="rowStart">The first row</param>
-        /// <param name="rowEnd">The last row</param>
+        /// <param name="row">The firts row</param>
         /// <param name="column">Column index</param>
         /// <returns>Generated vector</returns>
-        private static double[] GenerateColumn(Matrix<double> a, int rowStart, int rowEnd, int column)
+        private static double[] GenerateColumn(Matrix<double> a, int row, int column)
         {
-            var ru = rowEnd - rowStart + 1;
+            var ru = a.RowCount - row;
             var u = new double[ru];
 
-            for (var i = rowStart; i <= rowEnd; i++)
+            for (var i = row; i < a.RowCount; i++)
             {
-                u[i - rowStart] = a.At(i, rowStart);
-                a.At(i, rowStart, 0.0);
+                u[i - row] = a.At(i, row);
+                a.At(i, row, 0.0);
             }
 
             var norm = u.Sum(t => t * t);
             norm = Math.Sqrt(norm);
 
-            if (rowStart == rowEnd || norm == 0)
+            if (row == a.RowCount - 1 || norm == 0)
             {
-                a.At(rowStart, column, -u[0]);
+                a.At(row, column, -u[0]);
                 u[0] = Math.Sqrt(2.0);
                 return u;
             }
@@ -121,7 +121,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
                 scale *= -1.0;
             }
 
-            a.At(rowStart, column, -1.0 / scale);
+            a.At(row, column, -1.0 / scale);
 
             for (var i = 0; i < ru; i++)
             {
@@ -145,35 +145,42 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
         /// <param name="u">Work array</param>
         /// <param name="a">Q or R matrices</param>
         /// <param name="rowStart">The first row</param>
-        /// <param name="rowEnd">The last row</param>
+        /// <param name="rowDim">The last row</param>
         /// <param name="columnStart">The first column</param>
-        /// <param name="columnEnd">The last column</param>
-        private static void ComputeQR(double[] u, Matrix<double> a, int rowStart, int rowEnd, int columnStart, int columnEnd)
+        /// <param name="columnDim">The last column</param>
+        /// <param name="availableCores">Number of available CPUs</param>
+        private static void ComputeQR(double[] u, Matrix<double> a, int rowStart, int rowDim, int columnStart, int columnDim, int availableCores)
         {
-            if (rowEnd < rowStart || columnEnd < columnStart)
+            if (rowDim < rowStart || columnDim < columnStart)
             {
                 return;
             }
 
-            var v = new double[columnEnd - columnStart + 1];
-            for (var j = columnStart; j <= columnEnd; j++)
-            {
-                v[j - columnStart] = 0.0;
-            }
+            var tmpColCount = columnDim - columnStart;
 
-            for (var i = rowStart; i <= rowEnd; i++)
+            if ((availableCores > 1) && (tmpColCount > 200))
             {
-                for (var j = columnStart; j <= columnEnd; j++)
-                {
-                    v[j - columnStart] = v[j - columnStart] + (u[i - rowStart] * a.At(i, j));
-                }
-            }
+                var tmpSplit = columnStart + (tmpColCount / 2);
+                var tmpCores = availableCores / 2;
 
-            for (var i = rowStart; i <= rowEnd; i++)
+                CommonParallel.Invoke(
+                    () => ComputeQR(u, a, rowStart, rowDim, columnStart, tmpSplit, tmpCores),
+                    () => ComputeQR(u, a, rowStart, rowDim, tmpSplit, columnDim, tmpCores));
+            }
+            else
             {
-                for (var j = columnStart; j <= columnEnd; j++)
+                for (var j = columnStart; j < columnDim; j++)
                 {
-                    a.At(i, j, a.At(i, j) - (u[i - rowStart] * v[j - columnStart]));
+                    var scale = 0.0;
+                    for (var i = rowStart; i < rowDim; i++)
+                    {
+                        scale += u[i - rowStart] * a.At(i, j);
+                    }
+
+                    for (var i = rowStart; i < rowDim; i++)
+                    {
+                        a.At(i, j, a.At(i, j) - (u[i - rowStart] * scale));
+                    }
                 }
             }
         }
diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/UserCholesky.cs b/src/Numerics/LinearAlgebra/Single/Factorization/UserCholesky.cs
index 6035fa83..3bff9b46 100644
--- a/src/Numerics/LinearAlgebra/Single/Factorization/UserCholesky.cs
+++ b/src/Numerics/LinearAlgebra/Single/Factorization/UserCholesky.cs
@@ -33,6 +33,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
     using System;
     using Generic;
     using Properties;
+    using Threading;
 
     /// <summary>
     /// <para>A class which encapsulates the functionality of a Cholesky factorization for user matrices.</para>
@@ -67,32 +68,74 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
 
             // Create a new matrix for the Cholesky factor, then perform factorization (while overwriting).
             CholeskyFactor = matrix.Clone();
-            for (var j = 0; j < CholeskyFactor.RowCount; j++)
+            var tmpColumn = new float[CholeskyFactor.RowCount];
+
+            // Main loop - along the diagonal
+            for (var ij = 0; ij < CholeskyFactor.RowCount; ij++)
             {
-                var d = 0.0;
-                for (var k = 0; k < j; k++)
+                // "Pivot" element
+                var tmpVal = CholeskyFactor.At(ij, ij);
+
+                if (tmpVal > 0.0)
                 {
-                    var s = 0.0f;
-                    for (var i = 0; i < k; i++)
+                    tmpVal = (float)Math.Sqrt(tmpVal);
+                    CholeskyFactor.At(ij, ij, tmpVal);
+                    tmpColumn[ij] = tmpVal;
+
+                    // Calculate multipliers and copy to local column
+                    // Current column, below the diagonal
+                    for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
                     {
-                        s += CholeskyFactor.At(k, i) * CholeskyFactor.At(j, i);
+                        CholeskyFactor.At(i, ij, CholeskyFactor.At(i, ij) / tmpVal);
+                        tmpColumn[i] = CholeskyFactor.At(i, ij);
                     }
 
-                    s = (matrix.At(j, k) - s) / CholeskyFactor.At(k, k);
-                    CholeskyFactor.At(j, k, s);
-                    d += s * s;
+                    // Remaining columns, below the diagonal
+                    DoCholeskyStep(CholeskyFactor, CholeskyFactor.RowCount, ij + 1, CholeskyFactor.RowCount, tmpColumn, Control.NumberOfParallelWorkerThreads);
                 }
-
-                d = matrix.At(j, j) - d;
-                if (d <= 0.0)
+                else
                 {
                     throw new ArgumentException(Resources.ArgumentMatrixPositiveDefinite);
                 }
 
-                CholeskyFactor.At(j, j, (float)Math.Sqrt(d));
-                for (var k = j + 1; k < CholeskyFactor.RowCount; k++)
+                for (var i = ij + 1; i < CholeskyFactor.RowCount; i++)
                 {
-                    CholeskyFactor.At(j, k, 0.0f);
+                    CholeskyFactor.At(ij, i, 0.0f);
+                }
+            }
+        }
+
+        /// <summary>
+        /// Calculate Cholesky step
+        /// </summary>
+        /// <param name="data">Factor matrix</param>
+        /// <param name="rowDim">Number of rows</param>
+        /// <param name="firstCol">Column start</param>
+        /// <param name="colLimit">Total columns</param>
+        /// <param name="multipliers">Multipliers calculated previously</param>
+        /// <param name="availableCores">Number of available processors</param>
+        private static void DoCholeskyStep(Matrix<float> data, int rowDim, int firstCol, int colLimit, float[] multipliers, int availableCores)
+        {
+            var tmpColCount = colLimit - firstCol;
+
+            if ((availableCores > 1) && (tmpColCount > 200))
+            {
+                var tmpSplit = firstCol + (tmpColCount / 3);
+                var tmpCores = availableCores / 2;
+
+                CommonParallel.Invoke(
+                    () => DoCholeskyStep(data, rowDim, firstCol, tmpSplit, multipliers, tmpCores),
+                    () => DoCholeskyStep(data, rowDim, tmpSplit, colLimit, multipliers, tmpCores));
+            }
+            else
+            {
+                for (var j = firstCol; j < colLimit; j++)
+                {
+                    var tmpVal = multipliers[j];
+                    for (var i = j; i < rowDim; i++)
+                    {
+                        data.At(i, j, data.At(i, j) - (multipliers[i] * tmpVal));
+                    }
                 }
             }
         }
diff --git a/src/Numerics/LinearAlgebra/Single/Factorization/UserQR.cs b/src/Numerics/LinearAlgebra/Single/Factorization/UserQR.cs
index 72a1a947..11fe1b93 100644
--- a/src/Numerics/LinearAlgebra/Single/Factorization/UserQR.cs
+++ b/src/Numerics/LinearAlgebra/Single/Factorization/UserQR.cs
@@ -34,6 +34,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
     using System.Linq;
     using Generic;
     using Properties;
+    using Threading;
 
     /// <summary>
     /// <para>A class which encapsulates the functionality of the QR decomposition.</para>
@@ -76,13 +77,13 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
             var u = new float[minmn][];
             for (var i = 0; i < minmn; i++)
             {
-                u[i] = GenerateColumn(MatrixR, i, matrix.RowCount - 1, i);
-                ComputeQR(u[i], MatrixR, i, matrix.RowCount - 1, i + 1, matrix.ColumnCount - 1);
+                u[i] = GenerateColumn(MatrixR, i, i);
+                ComputeQR(u[i], MatrixR, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
             }
 
             for (var i = minmn - 1; i >= 0; i--)
             {
-                ComputeQR(u[i], MatrixQ, i, matrix.RowCount - 1, i, matrix.RowCount - 1);
+                ComputeQR(u[i], MatrixQ, i, matrix.RowCount, i, matrix.RowCount, Control.NumberOfParallelWorkerThreads);
             }
         }
 
@@ -90,27 +91,26 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
         /// Generate column from initial matrix to work array
         /// </summary>
         /// <param name="a">Initial matrix</param>
-        /// <param name="rowStart">The first row</param>
-        /// <param name="rowEnd">The last row</param>
+        /// <param name="row">The first row</param>
         /// <param name="column">Column index</param>
         /// <returns>Generated vector</returns>
-        private static float[] GenerateColumn(Matrix<float> a, int rowStart, int rowEnd, int column)
+        private static float[] GenerateColumn(Matrix<float> a, int row, int column)
         {
-            var ru = rowEnd - rowStart + 1;
+            var ru = a.RowCount - row;
             var u = new float[ru];
 
-            for (var i = rowStart; i <= rowEnd; i++)
+            for (var i = row; i < a.RowCount; i++)
             {
-                u[i - rowStart] = a.At(i, rowStart);
-                a.At(i, rowStart, 0.0f);
+                u[i - row] = a.At(i, row);
+                a.At(i, row, 0.0f);
             }
 
             var norm = u.Sum(t => t * t);
             norm = (float)Math.Sqrt(norm);
 
-            if (rowStart == rowEnd || norm == 0)
+            if (row == a.RowCount - 1 || norm == 0)
             {
-                a.At(rowStart, column, -u[0]);
+                a.At(row, column, -u[0]);
                 u[0] = (float)Math.Sqrt(2.0);
                 return u;
             }
@@ -121,7 +121,7 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
                 scale *= -1.0f;
             }
 
-            a.At(rowStart, column, -1.0f / scale);
+            a.At(row, column, -1.0f / scale);
 
             for (var i = 0; i < ru; i++)
             {
@@ -145,35 +145,42 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
         /// <param name="u">Work array</param>
         /// <param name="a">Q or R matrices</param>
         /// <param name="rowStart">The first row</param>
-        /// <param name="rowEnd">The last row</param>
+        /// <param name="rowDim">The last row</param>
         /// <param name="columnStart">The first column</param>
-        /// <param name="columnEnd">The last column</param>
-        private static void ComputeQR(float[] u, Matrix<float> a, int rowStart, int rowEnd, int columnStart, int columnEnd)
+        /// <param name="columnDim">The last column</param>
+        /// <param name="availableCores">Number of available CPUs</param>
+        private static void ComputeQR(float[] u, Matrix<float> a, int rowStart, int rowDim, int columnStart, int columnDim, int availableCores)
         {
-            if (rowEnd < rowStart || columnEnd < columnStart)
+            if (rowDim < rowStart || columnDim < columnStart)
             {
                 return;
             }
 
-            var v = new float[columnEnd - columnStart + 1];
-            for (var j = columnStart; j <= columnEnd; j++)
-            {
-                v[j - columnStart] = 0.0f;
-            }
+            var tmpColCount = columnDim - columnStart;
 
-            for (var i = rowStart; i <= rowEnd; i++)
+            if ((availableCores > 1) && (tmpColCount > 200))
             {
-                for (var j = columnStart; j <= columnEnd; j++)
-                {
-                    v[j - columnStart] = v[j - columnStart] + (u[i - rowStart] * a.At(i, j));
-                }
-            }
+                var tmpSplit = columnStart + (tmpColCount / 2);
+                var tmpCores = availableCores / 2;
 
-            for (var i = rowStart; i <= rowEnd; i++)
+                CommonParallel.Invoke(
+                    () => ComputeQR(u, a, rowStart, rowDim, columnStart, tmpSplit, tmpCores),
+                    () => ComputeQR(u, a, rowStart, rowDim, tmpSplit, columnDim, tmpCores));
+            }
+            else
             {
-                for (var j = columnStart; j <= columnEnd; j++)
+                for (var j = columnStart; j < columnDim; j++)
                 {
-                    a.At(i, j, a.At(i, j) - (u[i - rowStart] * v[j - columnStart]));
+                    var scale = 0.0f;
+                    for (var i = rowStart; i < rowDim; i++)
+                    {
+                        scale += u[i - rowStart] * a.At(i, j);
+                    }
+
+                    for (var i = rowStart; i < rowDim; i++)
+                    {
+                        a.At(i, j, a.At(i, j) - (u[i - rowStart] * scale));
+                    }
                 }
             }
         }