native: added QR factor

src/NativeWrappers/MKL/lapack.cpp

@@ -523,6 +523,7 @@ extern "C" {
 				if (i > j)
 				{
 					q[j * m + i] = r[j * m + i];
+					r[j * m + i] = 0.0f;
 				}
 			}
 		}
@@ -552,6 +553,7 @@ extern "C" {
 				if (i > j)
 				{
 					q[j * m + i] = r[j * m + i];
+					r[j * m + i] = 0.0;
 				}
 			}
 		}
@@ -573,6 +575,10 @@ extern "C" {
 	{
 		int info = 0;
 		CGEQRF(&m, &n, r, &m, tau, work, &len, &info);
+		
+		MKL_Complex8 zero;
+		zero.real = 0.0f;
+		zero.imag = 0.0f;
 
 		for (int i = 0; i < m; ++i)
 		{
@@ -581,6 +587,7 @@ extern "C" {
 				if (i > j)
 				{
 					q[j * m + i] = r[j * m + i];
+					r[j * m + i] = zero;
 				}
 			}
 		}
@@ -603,6 +610,10 @@ extern "C" {
 		int info = 0;
 		ZGEQRF(&m, &n, r, &m, tau, work, &len, &info);
 
+		MKL_Complex16 zero;
+		zero.real = 0.0;
+		zero.imag = 0.0;
+
 		for (int i = 0; i < m; ++i)
 		{
 			for (int j = 0; j < m && j < n; ++j)
@@ -610,6 +621,7 @@ extern "C" {
 				if (i > j)
 				{
 					q[j * m + i] = r[j * m + i];
+					r[j * m + i] = zero;
 				}
 			}
 		}

src/Numerics/Algorithms/LinearAlgebra/ILinearAlgebraProviderOfT.cs

@@ -316,13 +316,15 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// Computes the QR factorization of A.
         /// </summary>
         /// <param name="r">On entry, it is the M by N A matrix to factor. On exit,
-        /// it is overwritten with the R matrix of the QR factorization. </param>
+        /// it is overwritten with the R matrix of the QR factorization.</param>
         /// <param name="rowsR">The number of rows in the A matrix.</param>
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
         /// QR factorization.</param>
+        /// <param name="tau">A min(m,n) vector. On exit, contains additional information
+        /// to be used by the QR solve routine.</param>
         /// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
-        void QRFactor(T[] r, int rowsR, int columnsR, T[] q);
+        void QRFactor(T[] r, int rowsR, int columnsR, T[] q, T[] tau);
 
         /// <summary>
         /// Computes the QR factorization of A.
@@ -333,11 +335,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="q">On exit, A M by M matrix that holds the Q matrix of the 
         /// QR factorization.</param>
+        /// <param name="tau">A min(m,n) vector. On exit, contains additional information
+        /// to be used by the QR solve routine.</param>
         /// <param name="work">The work array. The array must have a length of at least N,
         /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
         /// work size value.</param>
         /// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
-        void QRFactor(T[] r, int rowsR, int columnsR, T[] q, T[] work);
+        void QRFactor(T[] r, int rowsR, int columnsR, T[] q, T[] tau, T[] work);
 
         /// <summary>
         /// Solves A*X=B for X using QR factorization of A.

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex.cs

@@ -1359,8 +1359,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
         /// QR factorization.</param>
+        /// <param name="tau">A min(m,n) vector. On exit, contains additional information
+        /// to be used by the QR solve routine.</param>
         /// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
-        public virtual void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q)
+        public virtual void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau)
         {
             if (r == null)
             {
@@ -1374,16 +1376,21 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (r.Length != rowsR * columnsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "r");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
+            }
+
+            if (tau.Length < Math.Min(rowsR, columnsR))
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
             }
 
             if (q.Length != rowsR * rowsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "q");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
             }
 
             var work = new Complex[rowsR * rowsR];
-            QRFactor(r, rowsR, columnsR, q, work);
+            QRFactor(r, rowsR, columnsR, q, tau, work);
         }
 
         /// <summary>
@@ -1395,11 +1402,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="q">On exit, A M by M matrix that holds the Q matrix of the 
         /// QR factorization.</param>
+        /// <param name="tau">A min(m,n) vector. On exit, contains additional information
+        /// to be used by the QR solve routine.</param>
         /// <param name="work">The work array. The array must have a length of at least N,
         /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
         /// work size value.</param>
         /// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
-        public virtual void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] work)
+        public virtual void QRFactor(Complex[] r, int rowsR, int columnsR, Complex[] q, Complex[] tau, Complex[] work)
         {
             if (r == null)
             {
@@ -1418,12 +1427,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (r.Length != rowsR * columnsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "r");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
+            }
+
+            if (tau.Length < Math.Min(rowsR, columnsR))
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
             }
 
             if (q.Length != rowsR * rowsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "q");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
             }
 
             if (work.Length < rowsR * rowsR)
@@ -1681,8 +1695,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <summary>
         /// Solves A*X=B for X using a previously QR factored matrix.
         /// </summary>
-        /// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(Complex[],int,int,Complex[])"/>.</param>
-        /// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex[],int,int,Complex[])"/>. </param>
+        /// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(Complex[],int,int,Complex[],Complex[])"/>.</param>
+        /// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex[],int,int,Complex[],Complex[])"/>. </param>
         /// <param name="rowsR">The number of rows in the A matrix.</param>
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="b">The B matrix.</param>

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Complex32.cs

@@ -1359,8 +1359,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
         /// QR factorization.</param>
+        /// <param name="tau">A min(m,n) vector. On exit, contains additional information
+        /// to be used by the QR solve routine.</param>
         /// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
-        public virtual void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q)
+        public virtual void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau)
         {
             if (r == null)
             {
@@ -1374,16 +1376,21 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (r.Length != rowsR * columnsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "r");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
+            }
+
+            if (tau.Length < Math.Min(rowsR, columnsR))
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
             }
 
             if (q.Length != rowsR * rowsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "q");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
             }
 
             var work = new Complex32[rowsR * rowsR];
-            QRFactor(r, rowsR, columnsR, q, work);
+            QRFactor(r, rowsR, columnsR, q, tau, work);
         }
 
         /// <summary>
@@ -1395,11 +1402,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="q">On exit, A M by M matrix that holds the Q matrix of the 
         /// QR factorization.</param>
+        /// <param name="tau">A min(m,n) vector. On exit, contains additional information
+        /// to be used by the QR solve routine.</param>
         /// <param name="work">The work array. The array must have a length of at least N,
         /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
         /// work size value.</param>
         /// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
-        public virtual void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] work)
+        public virtual void QRFactor(Complex32[] r, int rowsR, int columnsR, Complex32[] q, Complex32[] tau, Complex32[] work)
         {
             if (r == null)
             {
@@ -1418,12 +1427,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (r.Length != rowsR * columnsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "r");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
+            }
+
+            if (tau.Length < Math.Min(rowsR, columnsR))
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
             }
 
             if (q.Length != rowsR * rowsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "q");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
             }
 
             if (work.Length < rowsR * rowsR)
@@ -1681,8 +1695,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <summary>
         /// Solves A*X=B for X using a previously QR factored matrix.
         /// </summary>
-        /// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[])"/>.</param>
-        /// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[])"/>. </param>
+        /// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[],Complex32[])"/>.</param>
+        /// <param name="r">The R matrix obtained by calling <see cref="QRFactor(Complex32[],int,int,Complex32[],Complex32[])"/>. </param>
         /// <param name="rowsR">The number of rows in the A matrix.</param>
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="b">The B matrix.</param>

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Double.cs

@@ -1353,8 +1353,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
         /// QR factorization.</param>
+        /// <param name="tau">A min(m,n) vector. On exit, contains additional information
+        /// to be used by the QR solve routine.</param>
         /// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
-        public virtual void QRFactor(double[] r, int rowsR, int columnsR, double[] q)
+        public virtual void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau)
         {
             if (r == null)
             {
@@ -1368,16 +1370,21 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (r.Length != rowsR * columnsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "r");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
+            }
+
+            if (tau.Length < Math.Min(rowsR, columnsR))
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
             }
 
             if (q.Length != rowsR * rowsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "q");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
             }
 
             var work = new double[rowsR * rowsR];
-            QRFactor(r, rowsR, columnsR, q, work);
+            QRFactor(r, rowsR, columnsR, q, tau, work);
         }
 
         /// <summary>
@@ -1389,11 +1396,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="q">On exit, A M by M matrix that holds the Q matrix of the 
         /// QR factorization.</param>
+        /// <param name="tau">A min(m,n) vector. On exit, contains additional information
+        /// to be used by the QR solve routine.</param>
         /// <param name="work">The work array. The array must have a length of at least N,
         /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
         /// work size value.</param>
         /// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
-        public virtual void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] work)
+        public virtual void QRFactor(double[] r, int rowsR, int columnsR, double[] q, double[] tau, double[] work)
         {
             if (r == null)
             {
@@ -1412,12 +1421,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (r.Length != rowsR * columnsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "r");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
+            }
+
+            if (tau.Length < Math.Min(rowsR, columnsR))
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
             }
 
             if (q.Length != rowsR * rowsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "q");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
             }
 
             if (work.Length < rowsR * rowsR)
@@ -1676,8 +1690,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <summary>
         /// Solves A*X=B for X using a previously QR factored matrix.
         /// </summary>
-        /// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(double[],int,int,double[])"/>.</param>
-        /// <param name="r">The R matrix obtained by calling <see cref="QRFactor(double[],int,int,double[])"/>. </param>
+        /// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(double[],int,int,double[],double[])"/>.</param>
+        /// <param name="r">The R matrix obtained by calling <see cref="QRFactor(double[],int,int,double[],double[])"/>. </param>
         /// <param name="rowsR">The number of rows in the A matrix.</param>
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="b">The B matrix.</param>

src/Numerics/Algorithms/LinearAlgebra/ManagedLinearAlgebraProvider.Single.cs

@@ -1354,8 +1354,10 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
         /// QR factorization.</param>
+        /// <param name="tau">A min(m,n) vector. On exit, contains additional information
+        /// to be used by the QR solve routine.</param>
         /// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
-        public virtual void QRFactor(float[] r, int rowsR, int columnsR, float[] q)
+        public virtual void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau)
         {
             if (r == null)
             {
@@ -1369,16 +1371,21 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (r.Length != rowsR * columnsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "r");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
+            }
+
+            if (tau.Length < Math.Min(rowsR, columnsR))
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
             }
 
             if (q.Length != rowsR * rowsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "q");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
             }
 
             var work = new float[rowsR * rowsR];
-            QRFactor(r, rowsR, columnsR, q, work);
+            QRFactor(r, rowsR, columnsR, q, tau, work);
         }
 
         /// <summary>
@@ -1390,11 +1397,13 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="q">On exit, A M by M matrix that holds the Q matrix of the 
         /// QR factorization.</param>
+        /// <param name="tau">A min(m,n) vector. On exit, contains additional information
+        /// to be used by the QR solve routine.</param>
         /// <param name="work">The work array. The array must have a length of at least N,
         /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
         /// work size value.</param>
         /// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
-        public virtual void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] work)
+        public virtual void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau, float[] work)
         {
             if (r == null)
             {
@@ -1413,12 +1422,17 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
 
             if (r.Length != rowsR * columnsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "r");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
+            }
+
+            if (tau.Length < Math.Min(rowsR, columnsR))
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
             }
 
             if (q.Length != rowsR * rowsR)
             {
-                throw new ArgumentException(Resources.ArgumentArraysSameLength, "q");
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
             }
 
             if (work.Length < rowsR * rowsR)
@@ -1677,8 +1691,8 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra
         /// <summary>
         /// Solves A*X=B for X using a previously QR factored matrix.
         /// </summary>
-        /// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(float[],int,int,float[])"/>.</param>
-        /// <param name="r">The R matrix obtained by calling <see cref="QRFactor(float[],int,int,float[])"/>. </param>
+        /// <param name="q">The Q matrix obtained by calling <see cref="QRFactor(float[],int,int,float[],float[])"/>.</param>
+        /// <param name="r">The R matrix obtained by calling <see cref="QRFactor(float[],int,int,float[],float[])"/>. </param>
         /// <param name="rowsR">The number of rows in the A matrix.</param>
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="b">The B matrix.</param>

src/Numerics/Algorithms/LinearAlgebra/native.generic.include

@@ -517,11 +517,39 @@
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="q">On exit, A M by M matrix that holds the Q matrix of the
         /// QR factorization.</param>
+        /// <param name="tau">A min(m,n) vector. On exit, contains additional information
+        /// to be used by the QR solve routine.</param>
         /// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
         [SecuritySafeCritical]
-        public override void QRFactor(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q)
+        public override void QRFactor(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] tau)
         {
-            throw new NotImplementedException();
+            if (r == null)
+            {
+                throw new ArgumentNullException("r");
+            }
+
+            if (q == null)
+            {
+                throw new ArgumentNullException("q");
+            }
+
+            if (r.Length != rowsR * columnsR)
+            {
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
+            }
+
+            if (tau.Length < Math.Min(rowsR, columnsR))
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
+            }
+
+            if (q.Length != rowsR * rowsR)
+            {
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
+            }
+
+            var work = new <#=dataType#>[columnsR * Control.BlockSize];
+            SafeNativeMethods.<#=prefix#>_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length);
         }
 
         /// <summary>
@@ -533,14 +561,52 @@
         /// <param name="columnsR">The number of columns in the A matrix.</param>
         /// <param name="q">On exit, A M by M matrix that holds the Q matrix of the 
         /// QR factorization.</param>
+        /// <param name="tau">A min(m,n) vector. On exit, contains additional information
+        /// to be used by the QR solve routine.</param>
         /// <param name="work">The work array. The array must have a length of at least N,
         /// but should be N*blocksize. The blocksize is machine dependent. On exit, work[0] contains the optimal
         /// work size value.</param>
         /// <remarks>This is similar to the GEQRF and ORGQR LAPACK routines.</remarks>
         [SecuritySafeCritical]
-        public override void QRFactor(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] work)
+        public override void QRFactor(<#=dataType#>[] r, int rowsR, int columnsR, <#=dataType#>[] q, <#=dataType#>[] tau, <#=dataType#>[] work)
         {
-            throw new NotImplementedException();
+            if (r == null)
+            {
+                throw new ArgumentNullException("r");
+            }
+
+            if (q == null)
+            {
+                throw new ArgumentNullException("q");
+            }
+
+            if (work == null)
+            {
+                throw new ArgumentNullException("q");
+            }
+
+            if (r.Length != rowsR * columnsR)
+            {
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * columnsR"), "r");
+            }
+
+            if (tau.Length < Math.Min(rowsR, columnsR))
+            {
+                throw new ArgumentException(string.Format(Resources.ArrayTooSmall, "min(m,n)"), "tau");
+            }
+
+            if (q.Length != rowsR * rowsR)
+            {
+                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, "rowsR * rowsR"), "q");
+            }
+
+            if (work.Length < columnsR * Control.BlockSize)
+            {
+                work[0] = columnsR * Control.BlockSize;
+                throw new ArgumentException(Resources.WorkArrayTooSmall, "work");
+            }
+
+            SafeNativeMethods.<#=prefix#>_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length);
         }
 
         /// <summary>

src/Numerics/Algorithms/LinearAlgebra/safe.native.common.include

@@ -187,27 +187,39 @@ namespace MathNet.Numerics.Algorithms.LinearAlgebra.<#= namespaceSuffix #>
         internal static extern int z_lu_solve(int n, int nrhs, Complex[] a,  [In, Out] Complex[] b);
 
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
-		internal static extern int s_cholesky_solve(int n, int nrhs, float[] a, [In, Out] float[] b);
+        internal static extern int s_cholesky_solve(int n, int nrhs, float[] a, [In, Out] float[] b);
 
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
-		internal static extern int d_cholesky_solve(int n, int nrhs, double[] a, [In, Out] double[] b);
+        internal static extern int d_cholesky_solve(int n, int nrhs, double[] a, [In, Out] double[] b);
 
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
-		internal static extern int c_cholesky_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b);
+        internal static extern int c_cholesky_solve(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b);
 
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
-		internal static extern int z_cholesky_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b);
+        internal static extern int z_cholesky_solve(int n, int nrhs, Complex[] a, [In, Out] Complex[] b);
 
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
-		internal static extern int s_cholesky_solve_factored(int n, int nrhs, float[] a, float[] b);
+        internal static extern int s_cholesky_solve_factored(int n, int nrhs, float[] a, [In, Out] float[] b);
 
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
-		internal static extern int d_cholesky_solve_factored(int n, int nrhs, double[] a, double[] b);
+        internal static extern int d_cholesky_solve_factored(int n, int nrhs, double[] a, [In, Out] double[] b);
 
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
-		internal static extern int c_cholesky_solve_factored(int n, int nrhs, Complex32[] a, Complex32[] b);
+        internal static extern int c_cholesky_solve_factored(int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b);
 
         [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
-		internal static extern int z_cholesky_solve_factored(int n, int nrhs, Complex[] a, Complex[] b);
+        internal static extern int z_cholesky_solve_factored(int n, int nrhs, Complex[] a, [In, Out] Complex[] b);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int s_qr_factor(int m, int n, [In, Out] float[] r, [In, Out] float[] tau, [In, Out] float[] q, [In, Out] float[] work, int len);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int d_qr_factor(int m, int n, [In, Out] double[] r, [In, Out] double[] tau, [In, Out] double[] q, [In, Out] double[] work, int len);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int c_qr_factor(int m, int n, [In, Out] Complex32[] r, [In, Out] Complex32[] tau, [In, Out] Complex32[] q, [In, Out] Complex32[] work, int len);
+
+        [DllImport(DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)]
+        internal static extern int z_qr_factor(int m, int n, [In, Out] Complex[] r, [In, Out] Complex[] tau, [In, Out] Complex[] q, [In, Out] Complex[] work, int len);
 
         #endregion LAPACK

src/Numerics/Control.cs

@@ -46,6 +46,11 @@ namespace MathNet.Numerics
         private static int _numberOfThreads = Environment.ProcessorCount;
 #endif
 
+        /// <summary>
+        /// Initial block size for the native linear algebra provider.
+        /// </summary>
+        private static int _blockSize = 512;
+
         /// <summary>
         /// Initializes static members of the Control class.
         /// </summary>
@@ -60,7 +65,7 @@ namespace MathNet.Numerics
         /// <summary>
         /// Gets or sets a value indicating whether the distribution classes check validate each parameter.
         /// For the multivariate distributions this could involve an expensive matrix factorization.
-        /// The default setting of this property is true.
+        /// The default setting of this property is <c>true</c>.
         /// </summary>
         public static bool CheckDistributionParameters { get; set; }
 
@@ -109,5 +114,26 @@ namespace MathNet.Numerics
             }
 #endif
         }
+
+        /// <summary>
+        /// Gets or sets the the block size to use for the native linear
+        /// algebra provider.
+        /// </summary>
+        /// <value>The block size. Must be at least 32.</value>
+        public static int BlockSize
+        {
+            get
+            {
+                return _blockSize;
+            }
+
+            set
+            {
+                if (_blockSize > 31)
+                {
+                    _blockSize = value;
+                }
+            }
+        }
     }
 }

src/Numerics/LinearAlgebra/Complex/Factorization/DenseQR.cs

@@ -46,6 +46,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
     /// </remarks>
     public class DenseQR : QR
     {
+        /// <summary>
+        ///  Gets or sets Tau vector. Contains additional information on Q - used for native solver.
+        /// </summary>
+        public Complex[] Tau
+        {
+            get;
+            set;
+        }
+
         /// <summary>
         /// Initializes a new instance of the <see cref="DenseQR"/> class. This object will compute the
         /// QR factorization when the constructor is called and cache it's factorization.
@@ -67,7 +76,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex.Factorization
 
             MatrixR = matrix.Clone();
             MatrixQ = new DenseMatrix(matrix.RowCount);
-            Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data);
+            Tau = new Complex[Math.Min(matrix.RowCount, matrix.ColumnCount)];
+            Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data, Tau);
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Complex32/Factorization/DenseQR.cs

@@ -46,6 +46,15 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
     /// </remarks>
     public class DenseQR : QR
     {
+        /// <summary>
+        ///  Gets or sets Tau vector. Contains additional information on Q - used for native solver.
+        /// </summary>
+        public Complex32[] Tau
+        {
+            get;
+            set;
+        }
+
         /// <summary>
         /// Initializes a new instance of the <see cref="DenseQR"/> class. This object will compute the
         /// QR factorization when the constructor is called and cache it's factorization.
@@ -67,7 +76,8 @@ namespace MathNet.Numerics.LinearAlgebra.Complex32.Factorization
 
             MatrixR = matrix.Clone();
             MatrixQ = new DenseMatrix(matrix.RowCount);
-            Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data);
+            Tau = new Complex32[Math.Min(matrix.RowCount, matrix.ColumnCount)];
+            Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data, Tau);
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Double/Factorization/DenseQR.cs

@@ -45,6 +45,15 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
     /// </remarks>
     public class DenseQR : QR
     {
+        /// <summary>
+        ///  Gets or sets Tau vector. Contains additional information on Q - used for native solver.
+        /// </summary>
+        public double[] Tau
+        {
+            get;
+            set;
+        }
+        
         /// <summary>
         /// Initializes a new instance of the <see cref="DenseQR"/> class. This object will compute the
         /// QR factorization when the constructor is called and cache it's factorization.
@@ -54,11 +63,6 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
         /// <exception cref="ArgumentException">If <paramref name="matrix"/> row count is less then column count</exception>
         public DenseQR(DenseMatrix matrix)
         {
-            if (matrix == null)
-            {
-                throw new ArgumentNullException("matrix");
-            }
-
             if (matrix.RowCount < matrix.ColumnCount)
             {
                 throw new ArgumentException(Resources.ArgumentMatrixDimensions);
@@ -66,7 +70,8 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
 
             MatrixR = matrix.Clone();
             MatrixQ = new DenseMatrix(matrix.RowCount);
-            Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data);
+            Tau = new double[Math.Min(matrix.RowCount, matrix.ColumnCount)];
+            Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data, Tau);
         }
 
         /// <summary>

src/Numerics/LinearAlgebra/Double/Factorization/UserQR.cs

@@ -91,7 +91,7 @@ namespace MathNet.Numerics.LinearAlgebra.Double.Factorization
         /// Generate column from initial matrix to work array
         /// </summary>
         /// <param name="a">Initial matrix</param>
-        /// <param name="row">The firts row</param>
+        /// <param name="row">The first row</param>
         /// <param name="column">Column index</param>
         /// <returns>Generated vector</returns>
         private static double[] GenerateColumn(Matrix<double> a, int row, int column)

src/Numerics/LinearAlgebra/Single/Factorization/DenseQR.cs

@@ -45,6 +45,15 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
     /// </remarks>
     public class DenseQR : QR
     {
+        /// <summary>
+        ///  Gets or sets Tau vector. Contains additional information on Q - used for native solver.
+        /// </summary>
+        internal float[] Tau
+        {
+            get;
+            set;
+        }
+
         /// <summary>
         /// Initializes a new instance of the <see cref="DenseQR"/> class. This object will compute the
         /// QR factorization when the constructor is called and cache it's factorization.
@@ -66,7 +75,8 @@ namespace MathNet.Numerics.LinearAlgebra.Single.Factorization
 
             MatrixR = matrix.Clone();
             MatrixQ = new DenseMatrix(matrix.RowCount);
-            Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data);
+            Tau = new float[Math.Min(matrix.RowCount, matrix.ColumnCount)];
+            Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Data, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix)MatrixQ).Data, Tau);
         }
 
         /// <summary>

src/Numerics/Properties/Resources.Designer.cs

@@ -69,6 +69,15 @@ namespace MathNet.Numerics.Properties {
             }
         }
         
+        /// <summary>
+        ///   Looks up a localized string similar to The given array is the wrong length. Should be {0}..
+        /// </summary>
+        internal static string ArgumentArrayWrongLength {
+            get {
+                return ResourceManager.GetString("ArgumentArrayWrongLength", resourceCulture);
+            }
+        }
+        
         /// <summary>
         ///   Looks up a localized string similar to The argument must be between 0 and 1..
         /// </summary>

src/Numerics/Properties/Resources.resx

@@ -351,4 +351,7 @@
   <data name="ArrayTooSmall" xml:space="preserve">
     <value>The given array is too small. It must be at least {0} long.</value>
   </data>
+  <data name="ArgumentArrayWrongLength" xml:space="preserve">
+    <value>The given array is the wrong length. Should be {0}.</value>
+  </data>
 </root>

src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs

@@ -641,7 +641,152 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double
             AssertHelpers.AlmostEqual(b[5], 0, 14);
         }
 
+        [Test]
+        public void CanComputeQRFactorSquareMatrix()
+        {
+            var matrix = _matrices["Square3x3"];
+            var r = new double[matrix.RowCount * matrix.ColumnCount];
+            Array.Copy(matrix.Data, r, r.Length);
+
+            var tau = new double[3];
+            var q = new double[matrix.RowCount * matrix.RowCount];
+            Provider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau);
+
+            var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
+            var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r);
+            var a = mq * mr;
+
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var col = 0; col < matrix.ColumnCount; col++)
+                {
+                    AssertHelpers.AlmostEqual(matrix[row, col], a[row, col], 14);                    
+                }
+            }
+        }
+
+        [Test]
+        public void CanComputeQRFactorTallMatrix()
+        {
+            var matrix = _matrices["Tall3x2"];
+            var r = new double[matrix.RowCount * matrix.ColumnCount];
+            Array.Copy(matrix.Data, r, r.Length);
+
+            var tau = new double[3];
+            var q = new double[matrix.RowCount * matrix.RowCount];
+            Provider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau);
+
+            var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r);
+            var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
+            var a = mq * mr;
+
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var col = 0; col < matrix.ColumnCount; col++)
+                {
+                    AssertHelpers.AlmostEqual(matrix[row, col], a[row, col], 14);
+                }
+            }
+        }
+
+        [Test]
+        public void CanComputeQRFactorWideMatrix()
+        {
+            var matrix = _matrices["Wide2x3"];
+            var r = new double[matrix.RowCount * matrix.ColumnCount];
+            Array.Copy(matrix.Data, r, r.Length);
+
+            var tau = new double[3];
+            var q = new double[matrix.RowCount * matrix.RowCount];
+            Provider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau);
+
+            var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r);
+            var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
+            var a = mq * mr;
+
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var col = 0; col < matrix.ColumnCount; col++)
+                {
+                    AssertHelpers.AlmostEqual(matrix[row, col], a[row, col], 14);
+                }
+            }
+        }
 
+        [Test]
+        public void CanComputeQRFactorSquareMatrixWithWorkArray()
+        {
+            var matrix = _matrices["Square3x3"];
+            var r = new double[matrix.RowCount * matrix.ColumnCount];
+            Array.Copy(matrix.Data, r, r.Length);
+
+            var tau = new double[3];
+            var q = new double[matrix.RowCount * matrix.RowCount];
+            var work = new double[matrix.ColumnCount * Control.BlockSize];
+            Provider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
+
+            var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
+            var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r);
+            var a = mq * mr;
+
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var col = 0; col < matrix.ColumnCount; col++)
+                {
+                    AssertHelpers.AlmostEqual(matrix[row, col], a[row, col], 14);
+                }
+            }
+        }
+
+        [Test]
+        public void CanComputeQRFactorTallMatrixWithWorkArray()
+        {
+            var matrix = _matrices["Tall3x2"];
+            var r = new double[matrix.RowCount * matrix.ColumnCount];
+            Array.Copy(matrix.Data, r, r.Length);
+
+            var tau = new double[3];
+            var q = new double[matrix.RowCount * matrix.RowCount];
+            var work = new double[matrix.ColumnCount * Control.BlockSize];
+            Provider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
+
+            var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r);
+            var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
+            var a = mq * mr;
+
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var col = 0; col < matrix.ColumnCount; col++)
+                {
+                    AssertHelpers.AlmostEqual(matrix[row, col], a[row, col], 14);
+                }
+            }
+        }
+
+        [Test]
+        public void CanComputeQRFactorWideMatrixWithWorkArray()
+        {
+            var matrix = _matrices["Wide2x3"];
+            var r = new double[matrix.RowCount * matrix.ColumnCount];
+            Array.Copy(matrix.Data, r, r.Length);
+
+            var tau = new double[3];
+            var q = new double[matrix.RowCount * matrix.RowCount];
+            var work = new double[matrix.ColumnCount * Control.BlockSize];
+            Provider.QRFactor(r, matrix.RowCount, matrix.ColumnCount, q, tau, work);
+
+            var mr = new DenseMatrix(matrix.RowCount, matrix.ColumnCount, r);
+            var mq = new DenseMatrix(matrix.RowCount, matrix.RowCount, q);
+            var a = mq * mr;
+
+            for (var row = 0; row < matrix.RowCount; row++)
+            {
+                for (var col = 0; col < matrix.ColumnCount; col++)
+                {
+                    AssertHelpers.AlmostEqual(matrix[row, col], a[row, col], 14);
+                }
+            }
+        }
 
         /// <summary>
         /// Checks to see if a matrix and array contain the same values.