/// <summary>
        /// Initializes a new instance of the <see cref="UserQR"/> class. This object will compute the
        /// QR factorization when the constructor is called and cache it's factorization.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <param name="method">The QR factorization method to use.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        public static UserQR Create(Matrix<Complex> matrix, QRMethod method = QRMethod.Full)
        {
            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw Matrix.DimensionsDontMatch<ArgumentException>(matrix);
            }

            Matrix<Complex> q;
            Matrix<Complex> r;

            var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount);
            var u = new Complex[minmn][];

            if (method == QRMethod.Full)
            {
                r = matrix.Clone();
                q = Matrix<Complex>.Build.SameAs(matrix, matrix.RowCount, matrix.RowCount);

                for (var i = 0; i < matrix.RowCount; i++)
                {
                    q.At(i, i, 1.0f);
                }

                for (var i = 0; i < minmn; i++)
                {
                    u[i] = GenerateColumn(r, i, i);
                    ComputeQR(u[i], r, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.MaxDegreeOfParallelism);
                }

                for (var i = minmn - 1; i >= 0; i--)
                {
                    ComputeQR(u[i], q, i, matrix.RowCount, i, matrix.RowCount, Control.MaxDegreeOfParallelism);
                }
            }
            else
            {
                q = matrix.Clone();

                for (var i = 0; i < minmn; i++)
                {
                    u[i] = GenerateColumn(q, i, i);
                    ComputeQR(u[i], q, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.MaxDegreeOfParallelism);
                }

                r = q.SubMatrix(0, matrix.ColumnCount, 0, matrix.ColumnCount);
                q.Clear();

                for (var i = 0; i < matrix.ColumnCount; i++)
                {
                    q.At(i, i, 1.0f);
                }

                for (var i = minmn - 1; i >= 0; i--)
                {
                    ComputeQR(u[i], q, i, matrix.RowCount, i, matrix.ColumnCount, Control.MaxDegreeOfParallelism);
                }
            }

            return new UserQR(q, r, method);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="UserQR"/> class. This object will compute the
        /// QR factorization when the constructor is called and cache it's factorization.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <param name="method">The QR factorization method to use.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        public static UserQR Create(Matrix<double> matrix, QRMethod method = QRMethod.Full)
        {
            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw Matrix.DimensionsDontMatch<ArgumentException>(matrix);
            }

            Matrix<double> q;
            Matrix<double> r;

            var minmn = Math.Min(matrix.RowCount, matrix.ColumnCount);
            var u = new double[minmn][];

            if (method == QRMethod.Full)
            {
                r = matrix.Clone();
                q = matrix.CreateMatrix(matrix.RowCount, matrix.RowCount);

                for (var i = 0; i < matrix.RowCount; i++)
                {
                    q.At(i, i, 1.0);
                }

                for (var i = 0; i < minmn; i++)
                {
                    u[i] = GenerateColumn(r, i, i);
                    ComputeQR(u[i], r, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
                }

                for (var i = minmn - 1; i >= 0; i--)
                {
                    ComputeQR(u[i], q, i, matrix.RowCount, i, matrix.RowCount, Control.NumberOfParallelWorkerThreads);
                }
            }
            else
            {
                q = matrix.Clone();

                for (var i = 0; i < minmn; i++)
                {
                    u[i] = GenerateColumn(q, i, i);
                    ComputeQR(u[i], q, i, matrix.RowCount, i + 1, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
                }

                r = q.SubMatrix(0, matrix.ColumnCount, 0, matrix.ColumnCount);
                q.Clear();

                for (var i = 0; i < matrix.ColumnCount; i++)
                {
                    q.At(i, i, 1.0);
                }

                for (var i = minmn - 1; i >= 0; i--)
                {
                    ComputeQR(u[i], q, i, matrix.RowCount, i, matrix.ColumnCount, Control.NumberOfParallelWorkerThreads);
                }
            }

            return new UserQR(q, r, method);
        }

        /// <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.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <param name="method">The QR factorization method to use.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If <paramref name="matrix"/> row count is less then column count</exception>
        public static DenseQR Create(DenseMatrix matrix, QRMethod method = QRMethod.Full)
        {
            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw Matrix.DimensionsDontMatch<ArgumentException>(matrix);
            }

            var tau = new float[Math.Min(matrix.RowCount, matrix.ColumnCount)];
            Matrix<float> q;
            Matrix<float> r;

            if (method == QRMethod.Full)
            {
                r = matrix.Clone();
                q = new DenseMatrix(matrix.RowCount);
                Control.LinearAlgebraProvider.QRFactor(((DenseMatrix) r).Values, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix) q).Values, tau);
            }
            else
            {
                q = matrix.Clone();
                r = new DenseMatrix(matrix.ColumnCount);
                Control.LinearAlgebraProvider.ThinQRFactor(((DenseMatrix) q).Values, matrix.RowCount, matrix.ColumnCount, ((DenseMatrix) r).Values, tau);
            }

            return new DenseQR(q, r, method, tau);
        }

        /// <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.
        /// </summary>
        /// <param name="matrix">The matrix to factor.</param>
        /// <param name="method">The QR factorization method to use.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="matrix"/> is <c>null</c>.</exception>
        /// <exception cref="ArgumentException">If <paramref name="matrix"/> row count is less then column count</exception>
        public DenseQR(DenseMatrix matrix, QRMethod method = QRMethod.Full)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.RowCount < matrix.ColumnCount)
            {
                throw Matrix.DimensionsDontMatch<ArgumentException>(matrix);
            }

            Tau = new Complex32[Math.Min(matrix.RowCount, matrix.ColumnCount)];

            if (method == QRMethod.Full)
            {
                MatrixR = matrix.Clone();
                MatrixQ = new DenseMatrix(matrix.RowCount);
                Control.LinearAlgebraProvider.QRFactor(((DenseMatrix)MatrixR).Values, matrix.RowCount, matrix.ColumnCount,
                                                       ((DenseMatrix)MatrixQ).Values, Tau);
            }
            else
            {
                MatrixQ = matrix.Clone();
                MatrixR = new DenseMatrix(matrix.ColumnCount);
                Control.LinearAlgebraProvider.ThinQRFactor(((DenseMatrix)MatrixQ).Values, matrix.RowCount, matrix.ColumnCount,
                                                       ((DenseMatrix)MatrixR).Values, Tau);
            }
        }

 public override void QRSolveFactored(double[] q, double[] r, int rowsR, int columnsR, double[] tau, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full)
 {
     var work = new double[columnsR*Control.BlockSize];
     QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work, method);
 }

 public override void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full)
 {
     var work = new double[columns*Control.BlockSize];
     QRSolve(a, rows, columns, b, columnsB, x, work, method);
 }

        public void CanSolveForMatrixWithTallRandomMatrix(QRMethod method)
        {
            var matrixA = Matrix<Complex>.Build.Random(20, 10, 1);
            var matrixACopy = matrixA.Clone();
            var factorQR = matrixA.QR(method);

            var matrixB = Matrix<Complex>.Build.Random(20, 5, 1);
            var matrixX = factorQR.Solve(matrixB);

            // The solution X row dimension is equal to the column dimension of A
            Assert.AreEqual(matrixA.ColumnCount, matrixX.RowCount);

            // The solution X has the same number of columns as B
            Assert.AreEqual(matrixB.ColumnCount, matrixX.ColumnCount);

            var test = (matrixA.ConjugateTranspose() * matrixA).Inverse() * matrixA.ConjugateTranspose() * matrixB;

            for (var i = 0; i < matrixX.RowCount; i++)
            {
                for (var j = 0; j < matrixX.ColumnCount; j++)
                {
                    AssertHelpers.AlmostEqual(test[i, j], matrixX[i, j], 12);
                }
            }

            // Make sure A didn't change.
            for (var i = 0; i < matrixA.RowCount; i++)
            {
                for (var j = 0; j < matrixA.ColumnCount; j++)
                {
                    Assert.AreEqual(matrixACopy[i, j], matrixA[i, j]);
                }
            }
        }

        public override void QRFactor(float[] r, int rowsR, int columnsR, float[] q, float[] tau, QRMethod method = QRMethod.Full)
        {
            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 float[columnsR * Control.BlockSize];
            SafeNativeMethods.s_qr_factor(rowsR, columnsR, r, tau, q, work, work.Length);
        }

        public override void QRSolveFactored(double[] q, double[] r, int rowsR, int columnsR, double[] tau, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full)
        {
            if (r == null)
            {
                throw new ArgumentNullException("r");
            }

            if (q == null)
            {
                throw new ArgumentNullException("q");
            }

            if (b == null)
            {
                throw new ArgumentNullException("q");
            }

            if (x == null)
            {
                throw new ArgumentNullException("q");
            }

            if (r.Length != rowsR*columnsR)
            {
                throw new ArgumentException(Resources.ArgumentArraysSameLength, "r");
            }

            if (q.Length != rowsR*rowsR)
            {
                throw new ArgumentException(Resources.ArgumentArraysSameLength, "q");
            }

            if (b.Length != rowsR*columnsB)
            {
                throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
            }

            if (x.Length != columnsR*columnsB)
            {
                throw new ArgumentException(Resources.ArgumentArraysSameLength, "x");
            }

            if (rowsR < columnsR)
            {
                throw new ArgumentException(Resources.RowsLessThanColumns);
            }

            var work = new double[columnsR*Control.BlockSize];
            QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, work);
        }

        /// <summary>
        /// Solves A*X=B for X using QR factorization of A.
        /// </summary>
        /// <param name="a">The A matrix.</param>
        /// <param name="rows">The number of rows in the A matrix.</param>
        /// <param name="columns">The number of columns in the A matrix.</param>
        /// <param name="b">The B matrix.</param>
        /// <param name="columnsB">The number of columns of B.</param>
        /// <param name="x">On exit, the solution matrix.</param>
        /// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
        /// <remarks>Rows must be greater or equal to columns.</remarks>
        public virtual void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full)
        {
            if (a == null)
            {
                throw new ArgumentNullException("a");
            }

            if (b == null)
            {
                throw new ArgumentNullException("b");
            }

            if (x == null)
            {
                throw new ArgumentNullException("x");
            }

            if (a.Length != rows*columns)
            {
                throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
            }

            if (b.Length != rows*columnsB)
            {
                throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
            }

            if (x.Length != columns*columnsB)
            {
                throw new ArgumentException(Resources.ArgumentArraysSameLength, "x");
            }

            if (rows < columns)
            {
                throw new ArgumentException(Resources.RowsLessThanColumns);
            }

            var work = new float[rows * columns];

            var clone = new float[a.Length];
            a.Copy(clone);

            if (method == QRMethod.Full)
            {
                var q = new float[rows*rows];
                QRFactor(clone, rows, columns, q, work);
                QRSolveFactored(q, clone, rows, columns, null, b, columnsB, x, method);
            }
            else
            {
                var r = new float[columns*columns];
                ThinQRFactor(clone, rows, columns, r, work);
                QRSolveFactored(clone, r, rows, columns, null, b, columnsB, x, method);
            }
        }

 /// <summary>
 /// Solves A*X=B for X using a previously QR factored matrix.
 /// </summary>
 /// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
 /// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</param>
 /// <param name="r">The R matrix obtained by calling <see cref="QRFactor(double[],int,int,double[],double[])"/>. </param>
 /// <param name="rowsA">The number of rows in the A matrix.</param>
 /// <param name="columnsA">The number of columns in the A matrix.</param>
 /// <param name="tau">Contains additional information on Q. Only used for the native solver
 /// and can be <c>null</c> for the managed provider.</param>
 /// <param name="b">On entry the B matrix; on exit the X matrix.</param>
 /// <param name="columnsB">The number of columns of B.</param>
 /// <param name="x">On exit, the solution matrix.</param>
 /// <param name="work">The work array - only used in the native provider. 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>
 /// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
 /// <remarks>Rows must be greater or equal to columns.</remarks>
 public virtual void QRSolveFactored(double[] q, double[] r, int rowsA, int columnsA, double[] tau, double[] b, int columnsB, double[] x, double[] work, QRMethod method = QRMethod.Full)
 {
     QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, method);
 }

        public override void QRSolveFactored(float[] q, float[] r, int rowsA, int columnsA, float[] tau, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full)
        {
            if (r == null)
            {
                throw new ArgumentNullException("r");
            }

            if (q == null)
            {
                throw new ArgumentNullException("q");
            }

            if (b == null)
            {
                throw new ArgumentNullException("q");
            }

            if (x == null)
            {
                throw new ArgumentNullException("q");
            }

            int rowsQ, columnsQ, rowsR, columnsR;
            if (method == QRMethod.Full)
            {
                rowsQ = columnsQ = rowsR = rowsA;
                columnsR = columnsA;
            }
            else
            {
                rowsQ = rowsA;
                columnsQ = rowsR = columnsR = columnsA;
            }

            if (r.Length != rowsR * columnsR)
            {
                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsR * columnsR), "r");
            }

            if (q.Length != rowsQ * columnsQ)
            {
                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsQ * columnsQ), "q");
            }

            if (b.Length != rowsA * columnsB)
            {
                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, rowsA * columnsB), "b");
            }

            if (x.Length != columnsA * columnsB)
            {
                throw new ArgumentException(string.Format(Resources.ArgumentArrayWrongLength, columnsA * columnsB), "x");
            }

            if (method == QRMethod.Full)
            {
                var info = SafeNativeMethods.s_qr_solve_factored(rowsA, columnsA, columnsB, r, b, tau, x);
                
                if (info == (int)NativeError.MemoryAllocation)
                {
                    throw new MemoryAllocationException();
                }

                if (info < 0)
                {
                    throw new InvalidParameterException(Math.Abs(info));
                }
            }
            else
            {
                // we don't have access to the raw Q matrix any more(it is stored in R in the full QR), need to think about this.
                // let just call the managed version in the meantime. The heavy lifting has already been done. -marcus
                base.QRSolveFactored(q, r, rowsA, columnsA, tau, b, columnsB, x, QRMethod.Thin);
            }
        }

 /// <summary>
 /// Solves A*X=B for X using QR factorization of A.
 /// </summary>
 /// <param name="a">The A matrix.</param>
 /// <param name="rows">The number of rows in the A matrix.</param>
 /// <param name="columns">The number of columns in the A matrix.</param>
 /// <param name="b">The B matrix.</param>
 /// <param name="columnsB">The number of columns of B.</param>
 /// <param name="x">On exit, the solution matrix.</param>
 /// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
 /// <remarks>Rows must be greater or equal to columns.</remarks>
 public virtual void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full)
 {
     var work = new double[rows * columns];
     QRSolve(a, rows, columns, b, columnsB, x, work, method);
 }

 /// <summary>
 /// Solves A*X=B for X using a previously QR factored matrix.
 /// </summary>
 /// <param name="q">The Q matrix obtained by QR factor. This is only used for the managed provider and can be
 /// <c>null</c> for the native provider. The native provider uses the Q portion stored in the R matrix.</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="tau">Contains additional information on Q. Only used for the native solver
 /// and can be <c>null</c> for the managed provider.</param>
 /// <param name="b">On entry the B matrix; on exit the X matrix.</param>
 /// <param name="columnsB">The number of columns of B.</param>
 /// <param name="x">On exit, the solution matrix.</param>
 /// <param name="work">The work array - only used in the native provider. 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>
 /// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
 /// <remarks>Rows must be greater or equal to columns.</remarks>
 public virtual void QRSolveFactored(Complex[] q, Complex[] r, int rowsR, int columnsR, Complex[] tau, Complex[] b, int columnsB, Complex[] x, Complex[] work, QRMethod method = QRMethod.Full)
 {
     QRSolveFactored(q, r, rowsR, columnsR, tau, b, columnsB, x, method);
 }

 /// <summary>
 /// Solves A*X=B for X using QR factorization of A.
 /// </summary>
 /// <param name="a">The A matrix.</param>
 /// <param name="rows">The number of rows in the A matrix.</param>
 /// <param name="columns">The number of columns in the A matrix.</param>
 /// <param name="b">The B matrix.</param>
 /// <param name="columnsB">The number of columns of B.</param>
 /// <param name="x">On exit, the solution matrix.</param>
 /// <param name="method">The type of QR factorization to perform. <seealso cref="QRMethod"/></param>
 /// <remarks>Rows must be greater or equal to columns.</remarks>
 public virtual void QRSolve(Complex[] a, int rows, int columns, Complex[] b, int columnsB, Complex[] x, QRMethod method = QRMethod.Full)
 {
     var work = new Complex[rows*columns];
     QRSolve(a, rows, columns, b, columnsB, x, work, method);
 }

 public override QR<Complex32> QR(QRMethod method = QRMethod.Thin)
 {
     return DenseQR.Create(this, method);
 }

        public override void QRSolve(float[] a, int rows, int columns, float[] b, int columnsB, float[] x, QRMethod method = QRMethod.Full)
        {
            if (a == null)
            {
                throw new ArgumentNullException("a");
            }

            if (b == null)
            {
                throw new ArgumentNullException("b");
            }

            if (x == null)
            {
                throw new ArgumentNullException("x");
            }

            if (a.Length != rows * columns)
            {
                throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
            }

            if (b.Length != rows * columnsB)
            {
                throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
            }

            if (x.Length != columns * columnsB)
            {
                throw new ArgumentException(Resources.ArgumentArraysSameLength, "x");
            }

            if (rows < columns)
            {
                throw new ArgumentException(Resources.RowsLessThanColumns);
            }

            var info = SafeNativeMethods.s_qr_solve(rows, columns, columnsB, a, b, x);

            if (info == (int)NativeError.MemoryAllocation)
            {
                throw new MemoryAllocationException();
            }

            if (info < 0)
            {
                throw new InvalidParameterException(Math.Abs(info));
            }

            if (info > 0)
            {
                throw new ArgumentException(Resources.ArgumentMatrixNotRankDeficient, "a");
            }
        }

        /// <summary>
        /// Solves A*X=B for X using QR factorization of A.
        /// </summary>
        /// <param name="a">The A matrix.</param>
        /// <param name="rows">The number of rows in the A matrix.</param>
        /// <param name="columns">The number of columns in the A matrix.</param>
        /// <param name="b">The B matrix.</param>
        /// <param name="columnsB">The number of columns of B.</param>
        /// <param name="x">On exit, the solution matrix.</param>
        /// <remarks>Rows must be greater or equal to columns.</remarks>
        public override void QRSolve(double[] a, int rows, int columns, double[] b, int columnsB, double[] x, QRMethod method = QRMethod.Full)
        {
            if (a == null)
            {
                throw new ArgumentNullException("a");
            }

            if (b == null)
            {
                throw new ArgumentNullException("b");
            }

            if (x == null)
            {
                throw new ArgumentNullException("x");
            }

            if (a.Length != rows*columns)
            {
                throw new ArgumentException(Resources.ArgumentArraysSameLength, "a");
            }

            if (b.Length != rows*columnsB)
            {
                throw new ArgumentException(Resources.ArgumentArraysSameLength, "b");
            }

            if (x.Length != columns*columnsB)
            {
                throw new ArgumentException(Resources.ArgumentArraysSameLength, "x");
            }

            if (rows < columns)
            {
                throw new ArgumentException(Resources.RowsLessThanColumns);
            }

            var work = new double[columns*Control.BlockSize];
            QRSolve(a, rows, columns, b, columnsB, x, work);
        }

 DenseQR(Matrix<float> q, Matrix<float> rFull, QRMethod method, float[] tau)
     : base(q, rFull, method)
 {
     Tau = tau;
 }

 UserQR(Matrix<Complex> q, Matrix<Complex> rFull, QRMethod method)
     : base(q, rFull, method)
 {
 }