using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using UnitTestSharp;
using Azimuth.DenseLinearAlgebra;

namespace Azimuth.UnitTests.DenseLinearAlgebra
{
    public class QRDecompositionTests : TestFixture
    {
        public class ToHelpBuild : TestFixture
        {
            public void SingleColumnNiceNumbers()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 3, },
                    { 4, },
                });

                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    { .6, },
                    { .8, },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { 5.0, },
                });

                var expectedP = new int[] { 0 };

                var decomp = new QRDecomposition(matrix);

                var actualR = new RectangularMatrix(decomp.R);
                var actualQ = new RectangularMatrix(decomp.Q);

                CheckEqual(expectedQ, actualQ);
                CheckEqual(expectedR, actualR);
                CheckEqual(expectedP, decomp.P);
            }

            public void SingleColumnNiceNumbers_2()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    {  1, },
                    { -1, },
                });

                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    {  Math.Sqrt(2.0)/2.0, },
                    { -Math.Sqrt(2.0)/2.0, },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { Math.Sqrt(2.0), },
                });

                var expectedP = new int[] { 0 };

                var decomp = new QRDecomposition(matrix);

                var actualR = new RectangularMatrix(decomp.R);
                var actualQ = new RectangularMatrix(decomp.Q);

                CheckEqual(expectedQ, actualQ);
                CheckEqual(expectedR, actualR);
                CheckEqual(expectedP, decomp.P);
            }

            public void SingleRowNiceNumbers()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 3, 4, },
                });

                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    { 1, },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { 4, 3, },
                });

                var expectedP = new int[] { 1, 0 };

                var decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }
        }

        public class FullRankSquare : TestFixture
        {
            public void x11()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                var expectedQ = new SquareMatrix(new Scalar[,] {
                    { -1.0 },
                });

                var expectedR = new SquareMatrix(new Scalar[,] {
                    { 7.0 },
                });

                var expectedP = new int[] { 0 };

                QRDecomposition decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }

            public void x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {2, 1, },
                    {-1, 0, },
                });

                Scalar sqr5 = Math.Sqrt(5);
                
                var expectedQ = new SquareMatrix(new Scalar[,] {
                    {2/sqr5, 1/sqr5, },
                    {-1/sqr5, 2/sqr5, },
                });

                var expectedR = new SquareMatrix(new Scalar[,] {
                    {sqr5, 2*sqr5/5, },
                    {0,    sqr5/5, },
                });

                var expectedP = new int[] { 0, 1 };

                QRDecomposition decomp = new QRDecomposition(matrix.Clone());

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }

            public void x33()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { 2,  1,  3, },
                    {-1,  0,  7, },
                    { 0, -1, -1, },
                });

                Scalar sqr59 = Math.Sqrt(59);
                Scalar sqr17346 = Math.Sqrt(17346);
                Scalar sqr294 = Math.Sqrt(294);
                Scalar sqr4336p5 = Math.Sqrt(4336.5);
                
                var expectedQ = new SquareMatrix(new Scalar[,] {
                    {  3/sqr59,   121/sqr17346,  -1/sqr294, },
                    {  7/sqr59,   -52/sqr17346,  -2/sqr294, },
                    { -1/sqr59,    -1/sqr17346, -17/sqr294, },
                });

                var expectedR = new SquareMatrix(new Scalar[,] {
                    { sqr59, -Math.Sqrt(147.0/8673.0),              4/sqr59, },
                    {     0,             sqr294/sqr59,         61/sqr4336p5, },
                    {     0,                        0, Math.Sqrt(128.0/147), },
                });

                var expectedP = new int[] { 2, 0, 1, };

                QRDecomposition decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }

            public void Identity_x33()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { 1, 0, 0, },
                    { 0, 1, 0, },
                    { 0, 0, 1, },
                });

                var expectedQ = new SquareMatrix(new Scalar[,] {
                    { 1, 0, 0, },
                    { 0, 1, 0, },
                    { 0, 0, 1, },
                });

                var expectedR = new SquareMatrix(new Scalar[,] {
                    { 1, 0, 0, },
                    { 0, 1, 0, },
                    { 0, 0, 1, },
                });

                var expectedP = new int[] { 0, 1, 2, };

                QRDecomposition decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }
        }

        public class RankDeficientSquare : TestFixture
        {
            public void x22R1()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    {1, 2, },
                    {2, 4, },
                });

                Scalar sqrp2 = Math.Sqrt(.2);
                Scalar sqrp8 = Math.Sqrt(.8);
                Scalar sqr20 = Math.Sqrt(20);

                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    { sqrp2, },
                    { sqrp8, },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { sqr20, sqr20/2, },
                });

                var expectedP = new int[] { 1, 0, };

                QRDecomposition decomp = new QRDecomposition(matrix);
                
                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }

            public void x33R1()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {1, 2, 3, },
                    {2, 4, 6, },
                    {3, 6, 9, },
                });

                Scalar sqr14 = Math.Sqrt(14);
                Scalar sqr126 = Math.Sqrt(126);
                Scalar sqr56 = Math.Sqrt(56);

                //Only one column is determined
                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    {   sqr14/14, },
                    { 2*sqr14/14, },
                    { 3*sqr14/14, },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { sqr126, sqr14, sqr56, },
                });

                var expectedP = new int[] { 2, 0, 1 };

                QRDecomposition decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }

            public void x33R2()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { 13, 18, 48,  },
                    { 18, 25, 66,  },
                    { 48, 66, 180, },
                });

                Scalar sqr2797 = Math.Sqrt(2797);
                Scalar sqr181 = Math.Sqrt(181);
                Scalar sqr1085 = Math.Sqrt(1085);

                //Only one column is determined
                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    {              8/sqr1085,  169/sqr1085/sqr181, },
                    {             11/sqr1085,  368/sqr1085/sqr181, },
                    { Math.Sqrt(180.0/217.0), -180/sqr1085/sqr181, },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    {  6510/sqr1085,     2399/sqr1085,   1742/sqr1085, },
                    {             0, 2*sqr181/sqr1085, sqr181/sqr1085, },
                });

                var expectedP = new int[] { 2, 1, 0 };

                QRDecomposition decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }

            public void x33R2_0Cross()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {  9, 0,  7,  },
                    {  0, 0,  0,  },
                    {  7, 0, 25, },
                }); 

                Scalar sqr674 = Math.Sqrt(674);
                
                //Only one column is determined
                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    {  7/sqr674, 25/sqr674, },
                    {         0,         0, },
                    { 25/sqr674, -7/sqr674, },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { sqr674, 238/sqr674, 0, },
                    {      0, 176/sqr674, 0, },
                });

                var expectedP = new int[] { 2, 0, 1, };

                QRDecomposition decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }
        }

        public class FullRankTall : TestFixture
        {
            public void x42()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, },
                    { 2, 4, },
                    { 3, 6, },
                    { 4, 8, },
                });

                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    { 1.0/Math.Sqrt(30), },
                    { 2.0/Math.Sqrt(30), },
                    { 3.0/Math.Sqrt(30), },
                    { 4.0/Math.Sqrt(30), },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { 2.0*Math.Sqrt(30), 1.0*Math.Sqrt(30), },
                });

                var expectedP = new int[] { 1, 0, };

                QRDecomposition decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }

            public void x63()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 4, 1, },
                    { 2, 3, 0, },
                    { 3, 2, 0, },
                    { 4, 1, 1, },
                });

                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    { 1.0/Math.Sqrt(30),  2.0/Math.Sqrt(6),  0.5, },
                    { 2.0/Math.Sqrt(30),  1.0/Math.Sqrt(6), -0.5, },
                    { 3.0/Math.Sqrt(30),                 0, -0.5, },
                    { 4.0/Math.Sqrt(30), -1.0/Math.Sqrt(6),  0.5, },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { Math.Sqrt(30), 1.0 / Math.Sqrt(.075), 1.0/Math.Sqrt(1.2), },
                    {             0, 30.0 / Math.Sqrt(54),  1.0/Math.Sqrt(6), },
                    {             0,                   0,                  1, },
                });

                var expectedP = new int[] { 0, 1, 2, };

                QRDecomposition decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }
        }

        public class RankDeficientTall : TestFixture
        {
            public void x42()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, },
                    { 2, 4, },
                    { 3, 6, },
                    { 4, 8, },
                });

                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    { 1.0/Math.Sqrt(30), },
                    { 2.0/Math.Sqrt(30), },
                    { 3.0/Math.Sqrt(30), },
                    { 4.0/Math.Sqrt(30), },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { 2.0*Math.Sqrt(30.0), Math.Sqrt(30), },
                });

                var expectedP = new int[] { 1, 0, };

                QRDecomposition decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }

            public void x63()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, 3, },
                    { 2, 4, 6, },
                    { 3, 6, 9, },
                    { 4, 8, 12, },
                });

                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    { 1.0/Math.Sqrt(30), },
                    { 2.0/Math.Sqrt(30), },
                    { 3.0/Math.Sqrt(30), },
                    { 4.0/Math.Sqrt(30), },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { 3.0*Math.Sqrt(30.0), Math.Sqrt(30.0), 2.0*Math.Sqrt(30), },
                });

                var expectedP = new int[] { 2, 0, 1, };

                QRDecomposition decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }
        }

        public class RankDeficientFat : TestFixture
        {
            public void x24()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, 3, 4, },
                    { 2, 4, 6, 8, },
                });

                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    { 1.0/Math.Sqrt(5), },
                    { 2.0/Math.Sqrt(5), },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { 4.0*Math.Sqrt(5), Math.Sqrt(5), 2.0*Math.Sqrt(5), 3.0*Math.Sqrt(5), },
                });

                var expectedP = new int[] { 3, 0, 1, 2, };

                QRDecomposition decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }

            public void x36()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, 3, 4, },
                    { 2, 4, 6, 8, },
                    { 3, 6, 9, 12, },
                });

                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    { 1.0/Math.Sqrt(14), },
                    { 2.0/Math.Sqrt(14), },
                    { 3.0/Math.Sqrt(14), },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { 4.0*Math.Sqrt(14), Math.Sqrt(14), 2.0*Math.Sqrt(14), 3.0*Math.Sqrt(14), },
                });

                var expectedP = new int[] { 3, 0, 1, 2, };

                QRDecomposition decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }
        }

        public class FullRankFat : TestFixture
        {
            public void x24()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, 3, 4, },
                    { 2, 4, 6, 8, },
                });

                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    { 1.0/Math.Sqrt(5), },
                    { 2.0/Math.Sqrt(5), },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { 4.0*Math.Sqrt(5), Math.Sqrt(5), 2.0*Math.Sqrt(5), 3.0*Math.Sqrt(5), },
                });

                var expectedP = new int[] { 3, 0, 1, 2, };

                QRDecomposition decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }

            public void x36()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, 3, 4, },
                    { 4, 3, 2, 1, },
                    { 1, 0, 0, 1, },
                });

                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    { 1.0/Math.Sqrt(18),  7.0/Math.Sqrt(54),  1.0/Math.Sqrt(27), },
                    { 4.0/Math.Sqrt(18), -2.0/Math.Sqrt(54),  1.0/Math.Sqrt(27), },
                    { 1.0/Math.Sqrt(18),  1.0/Math.Sqrt(54), -5.0/Math.Sqrt(27), },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { Math.Sqrt(18),  9.0 / Math.Sqrt(18), 14.0/Math.Sqrt(18), 11.0 / Math.Sqrt(18), },
                    {             0, 27.0 / Math.Sqrt(54),  8.0/Math.Sqrt(54), 17.0 / Math.Sqrt(54), },
                    {             0,                   0,   5.0/Math.Sqrt(27),  5.0 / Math.Sqrt(27), },
                });

                var expectedP = new int[] { 0, 3, 1, 2, };

                QRDecomposition decomp = new QRDecomposition(matrix);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }
        }

        public class PostmultiplyByQ : TestFixture
        {
            public void WrongSize()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3, 4 });

                CheckThrow(typeof(Exception));
                decomp.PostmultiplyByQ(inputVector);
            }

            public void x11()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3 });
                var length = inputVector.Length;
                inputVector = decomp.PostmultiplyByQ(inputVector);

                CheckEqual(length, inputVector.Length);
                CheckEqual(3, inputVector[0]);
            }

            public void x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {2, 1, },
                    {-1, 0, },
                });

                Scalar sqr5 = Math.Sqrt(5);

                var expectedQ = new SquareMatrix(new Scalar[,] {
                    {2/sqr5, 1/sqr5, },
                    {-1/sqr5, 2/sqr5, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3, 4 });
                var outputVector1 = inputVector.Clone();
                var expectedVector = new DenseVector(new Scalar[] { -10.0 / sqr5, 5.0 / sqr5 });

                outputVector1 = decomp.PostmultiplyByQ(outputVector1);

                var outputVector2 = new DenseVector(decomp.Q.RowCount);
                outputVector2.AssignmentPostMatrixMultiply(inputVector, decomp.Q);
                
                CheckEqual(inputVector.Length, outputVector1.Length);
                CheckEqual(inputVector.Length, outputVector2.Length);

                CheckEqual(inputVector.Magnitude(), outputVector1.Magnitude());
                CheckEqual(inputVector.Magnitude(), outputVector2.Magnitude());

                CheckEqual(expectedVector, outputVector1);
                CheckEqual(expectedVector, outputVector2);
            }

            public void x42()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, },
                    { 2, 4, },
                    { 3, 6, },
                    { 4, 8, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { 1, 2, 3, 4 });
                var expected = new DenseVector(new Scalar[] { Math.Sqrt(30.0) });
                var outputVector = inputVector.Clone();

                outputVector = decomp.PostmultiplyByQ(outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x24()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, 3, 4, },
                    { 2, 4, 6, 8, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { Math.Sqrt(5.0), 2 * Math.Sqrt(5.0), });
                var expected = new DenseVector(new Scalar[] { 5.0 });
                var outputVector = inputVector.Clone();

                outputVector = decomp.PostmultiplyByQ(outputVector);

                CheckEqual(expected, outputVector);
            }
        }

        public class PremultiplyByQ : TestFixture
        {
            public void WrongSize()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3, 4 });

                CheckThrow(typeof(Exception));
                decomp.PremultiplyByQ(ref inputVector);
            }

            public void x11()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3 });
                var length = inputVector.Length;
                decomp.PremultiplyByQ(ref inputVector);

                CheckEqual(length, inputVector.Length);
                CheckEqual(3, inputVector[0]);
            }

            public void x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    {2, 1, },
                    {-1, 0, },
                });

                Scalar sqr5 = Math.Sqrt(5);

                var expectedQ = new SquareMatrix(new Scalar[,] {
                    {2/sqr5, 1/sqr5, },
                    {-1/sqr5, 2/sqr5, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3, 4 });
                var outputVector1 = inputVector.Clone();
                var expectedVector = new DenseVector(new Scalar[] { -2.0 / sqr5, 11.0 / sqr5 }); 

                decomp.PremultiplyByQ(ref outputVector1);

                var outputVector2 = new DenseVector(decomp.Q.ColumnCount);
                outputVector2.AssignmentPreMatrixMultiply(decomp.Q, inputVector);

                CheckEqual(inputVector.Length, outputVector1.Length);
                CheckEqual(inputVector.Length, outputVector2.Length);

                CheckEqual(inputVector.Magnitude(), outputVector1.Magnitude());
                CheckEqual(inputVector.Magnitude(), outputVector2.Magnitude());

                CheckEqual(expectedVector, outputVector1);
                CheckEqual(expectedVector, outputVector2);
            }

            public void x42()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, },
                    { 2, 4, },
                    { 3, 6, },
                    { 4, 8, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var expected = new DenseVector(new Scalar[] { 1, 2, 3, 4 });
                var inputVector = new DenseVector(new Scalar[] { Math.Sqrt(30.0) });
                var outputVector = new DenseVector(4);
                outputVector.SetLength(inputVector.Length);
                outputVector.Assignment(inputVector);

                decomp.PremultiplyByQ(ref outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x24()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, 3, 4, },
                    { 2, 4, 6, 8, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var expected = new DenseVector(new Scalar[] { Math.Sqrt(5.0), 2 * Math.Sqrt(5.0), });
                var inputVector = new DenseVector(new Scalar[] { 5.0 });
                var outputVector = new DenseVector(2);
                outputVector.SetLength(1);
                outputVector.Assignment(inputVector);

                decomp.PremultiplyByQ(ref outputVector);

                CheckEqual(expected, outputVector);
            }
        }

        public class PremultiplyByR : TestFixture
        {
            public void WrongSize()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3, 4 });

                CheckThrow(typeof(Exception));
                decomp.PremultiplyByR(inputVector);
            }

            public void x11()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3 });
                var length = inputVector.Length;
                inputVector = decomp.PremultiplyByR(inputVector);

                CheckEqual(length, inputVector.Length);
                CheckEqual(-21, inputVector[0]);
            }

            public void x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { 10, 1, },
                    { 0, 7, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var expected = new DenseVector(new Scalar[] { 12, 14, });
                var inputVector = new DenseVector(new Scalar[] { 1, 2 });
                var outputVector = inputVector.Clone();

                outputVector = decomp.PremultiplyByR(outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x42()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, },
                    { 2, 4, },
                    { 3, 6, },
                    { 4, 8, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { 1, 2, });

                var expected = new DenseVector(new Scalar[] { 4.0 * Math.Sqrt(30) });
                
                var outputVector = inputVector.Clone();

                outputVector = decomp.PremultiplyByR(outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x24()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, 3, 4, },
                    { 2, 4, 6, 8, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var hardCodedR = new DenseVector(new Scalar[] { 
                    Math.Sqrt(80.0),
                    Math.Sqrt(5.0),
                    Math.Sqrt(20.0),
                    Math.Sqrt(45.0),
                });
                var inputVector = new DenseVector(new Scalar[] { 1, 2, 3, 4, });
                var expected = new DenseVector(new Scalar[] { hardCodedR.DotProduct(inputVector) });
                var outputVector = inputVector.Clone();

                outputVector = decomp.PremultiplyByR(outputVector);

                CheckEqual(expected, outputVector);
            }
        }

        public class PostmultiplyByR : TestFixture
        {
            public void WrongSize()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3, 4 });

                CheckThrow(typeof(Exception));
                decomp.PostmultiplyByR(inputVector);
            }

            public void x11()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3 });
                var length = inputVector.Length;
                inputVector = decomp.PostmultiplyByR(inputVector);

                CheckEqual(length, inputVector.Length);
                CheckEqual(-21, inputVector[0]);
            }

            public void x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { 10, 1, },
                    { 0, 7, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var expected = new DenseVector(new Scalar[] { 10, 15, });
                var inputVector = new DenseVector(new Scalar[] { 1, 2 });
                var outputVector = inputVector.Clone();

                outputVector = decomp.PostmultiplyByR(outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x42()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, },
                    { 2, 4, },
                    { 3, 6, },
                    { 4, 8, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { Math.Sqrt(30) });

                var expected = new DenseVector(new Scalar[] { 60, 30 });

                var outputVector = inputVector.Clone();

                outputVector = decomp.PostmultiplyByR(outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x24()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, 3, 4, },
                    { 2, 4, 6, 8, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var hardCodedR = new DenseVector(new Scalar[] { 
                    Math.Sqrt(80.0),
                    Math.Sqrt(5.0),
                    Math.Sqrt(20.0),
                    Math.Sqrt(45.0),
                });
                var inputVector = new DenseVector(new Scalar[] { Math.Sqrt(5.0) });
                var expected = new DenseVector(new Scalar[] { 5 * 4, 5 * 1, 5 * 2, 5 * 3 });
                var outputVector = inputVector.Clone();

                outputVector = decomp.PostmultiplyByR(outputVector);

                CheckEqual(expected, outputVector);
            }
        }

        public class PremultiplyByInverseR_Tests : TestFixture
        {
            public void WrongSize()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3, 4 });

                CheckThrow(typeof(Exception));
                decomp.PremultiplyByInverseR(ref inputVector);
            }

            public void x11()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3 });
                var length = inputVector.Length;
                decomp.PremultiplyByInverseR(ref inputVector);

                CheckEqual(length, inputVector.Length);
                CheckEqual(-3.0/7.0, inputVector[0]);
            }

            public void x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { 10, 1, },
                    { 0, 7, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var expected = new DenseVector(new Scalar[] { 1.0/14.0, 4.0/14.0, });
                var inputVector = new DenseVector(new Scalar[] { 1, 2 });
                var outputVector = inputVector.Clone();

                decomp.PremultiplyByInverseR(ref outputVector);

                CheckEqual(expected, outputVector);
            }

            public void Pivoting()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { 1, 10, },
                    { 7, 0, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var expected = new DenseVector(new Scalar[] { 1.0 / 14.0, 4.0 / 14.0, });
                var inputVector = new DenseVector(new Scalar[] { 1, 2 });
                var outputVector = inputVector.Clone();

                decomp.PremultiplyByInverseR(ref outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x42()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 1, },
                    { 2, 1, },
                    { 3, 1, },
                    { 4, 1, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { 1, 2, });

                var expected = new DenseVector(new Scalar[] { -0.633922395092671, Math.Sqrt(6), });

                var outputVector = inputVector.Clone();

                decomp.PremultiplyByInverseR(ref outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x24()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, 3, 4, },
                    { 1, 1, 1, 1, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var hardCodedR = new DenseVector(new Scalar[] { 
                    Math.Sqrt(80.0),
                    Math.Sqrt(5.0),
                    Math.Sqrt(20.0),
                    Math.Sqrt(45.0),
                });
                var inputVector = new DenseVector(new Scalar[] { 1, 2, 3, 4, });
                var expected = new DenseVector(new Scalar[] { hardCodedR.DotProduct(inputVector) });
                var outputVector = inputVector.Clone();

                CheckThrow(typeof(Exception));
                decomp.PremultiplyByInverseR(ref outputVector);

                CheckEqual(expected, outputVector);
            }
        }

        public class PostmultiplyByInverseR_Tests : TestFixture
        {
            public void WrongSize()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3, 4 });

                CheckThrow(typeof(Exception));
                decomp.PostmultiplyByInverseR(ref inputVector);
            }

            public void x11()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3 });
                var length = inputVector.Length;
                decomp.PostmultiplyByInverseR(ref inputVector);

                CheckEqual(length, inputVector.Length);
                CheckEqual(-3.0/7.0, inputVector[0]);
            }

            public void x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { 10, 1, },
                    { 0, 7, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var expected = new DenseVector(new Scalar[] { 1.0/10.0, 19.0/70.0, });
                var inputVector = new DenseVector(new Scalar[] { 1, 2 });
                var outputVector = inputVector.Clone();

                decomp.PostmultiplyByInverseR(ref outputVector);

                CheckEqual(expected, outputVector);
            }

            public void Pivoting()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { 1, 10, },
                    { 7, 0, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var expected = new DenseVector(new Scalar[] { 1.0 / 10.0, 19.0 / 70.0, });
                var inputVector = new DenseVector(new Scalar[] { 1, 2 });
                var outputVector = inputVector.Clone();

                decomp.PostmultiplyByInverseR(ref outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x42()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 1, },
                    { 2, 1, },
                    { 3, 1, },
                    { 4, 1, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { 1, 2, });

                var expected = new DenseVector(new Scalar[] { 1.0/Math.Sqrt(30), 5.0/Math.Sqrt(6.0)});

                var outputVector = inputVector.Clone();

                decomp.PostmultiplyByInverseR(ref outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x24()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, 3, 4, },
                    { 1, 1, 1, 1, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var hardCodedR = new DenseVector(new Scalar[] { 
                    Math.Sqrt(80.0),
                    Math.Sqrt(5.0),
                    Math.Sqrt(20.0),
                    Math.Sqrt(45.0),
                });
                var inputVector = new DenseVector(new Scalar[] { Math.Sqrt(5.0) });
                var expected = new DenseVector(new Scalar[] { 5 * 4, 5 * 1, 5 * 2, 5 * 3 });
                var outputVector = inputVector.Clone();

                CheckThrow(typeof(Exception));
                decomp.PostmultiplyByInverseR(ref outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x33()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 3, 2, },
                    { 2, 4, 3, },
                    { 1, 5, 1, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { 11, 7, -3, });

                var expected = new DenseVector(new Scalar[] 
                {
                    1.55563491860281,
                    1.04903185661845,
                    -24.3563037459,
                });

                var outputVector = inputVector.Clone();

                decomp.PostmultiplyByInverseR(ref outputVector);

                CheckEqual(expected, outputVector);
            }
        }

        public class PostmultiplyByMatrix : TestFixture
        {
            public void WrongSize()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3, 4 });

                CheckThrow(typeof(Exception));
                decomp.PostmultiplyByMatrix(inputVector);
            }

            public void x11()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3 });
                var length = inputVector.Length;
                inputVector = decomp.PostmultiplyByMatrix(inputVector);

                CheckEqual(length, inputVector.Length);
                CheckEqual(21, inputVector[0]);
            }

            public void x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { 10, 1, },
                    { 0, 7, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var expected = new DenseVector(new Scalar[] { 10, 15, });
                var inputVector = new DenseVector(new Scalar[] { 1, 2 });
                var outputVector = inputVector.Clone();

                outputVector = decomp.PostmultiplyByMatrix(outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x42()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, },
                    { 2, 4, },
                    { 3, 6, },
                    { 4, 8, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { 1, 1, 1, 1, });

                var expected = new DenseVector(new Scalar[] { 10, 20 });

                var outputVector = inputVector.Clone();

                outputVector = decomp.PostmultiplyByMatrix(outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x24()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, 3, 4, },
                    { 2, 4, 6, 8, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { 2, 1 });
                var expected = new DenseVector(new Scalar[] { 4, 8, 12, 16 });
                var outputVector = inputVector.Clone();

                outputVector = decomp.PostmultiplyByMatrix(outputVector);

                CheckEqual(expected, outputVector);
            }
        }

        public class PremultiplyByMatrix : TestFixture
        {
            public void WrongSize()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3, 4 });

                CheckThrow(typeof(Exception));
                decomp.PremultiplyByMatrix(inputVector);
            }

            public void x11()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { -7 },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { -3 });
                var length = inputVector.Length;
                inputVector = decomp.PremultiplyByMatrix(inputVector);

                CheckEqual(length, inputVector.Length);
                CheckEqual(21, inputVector[0]);
            }

            public void x22()
            {
                var matrix = new SquareMatrix(new Scalar[,] {
                    { 10, 1, },
                    { 0, 7, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var expected = new DenseVector(new Scalar[] { 12, 14, });
                var inputVector = new DenseVector(new Scalar[] { 1, 2 });
                var outputVector = inputVector.Clone();

                outputVector = decomp.PremultiplyByMatrix(outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x42()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, },
                    { 2, 4, },
                    { 3, 6, },
                    { 4, 8, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { 1, 1, });

                var expected = new DenseVector(new Scalar[] { 3, 6, 9, 12 });

                var outputVector = inputVector.Clone();

                outputVector = decomp.PremultiplyByMatrix(outputVector);

                CheckEqual(expected, outputVector);
            }

            public void x24()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 2, 3, 4, },
                    { 2, 4, 6, 8, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);

                var inputVector = new DenseVector(new Scalar[] { 1, 1, 1, 1, });
                var expected = new DenseVector(new Scalar[] { 10, 20, });
                var outputVector = inputVector.Clone();

                outputVector = decomp.PremultiplyByMatrix(outputVector);

                CheckEqual(expected, outputVector);
            }
        }

        public class BasicSolution : TestFixture
        {
            public class FullRankSquare : TestFixture
            {
                public void x11()
                {
                    var matrix = new SquareMatrix(new Scalar[,] {
                        { -7 },
                    });

                    var inputOutput = new DenseVector(new Scalar[] {
                        7,
                    });

                    var expected = new DenseVector(new Scalar[] {
                        -1,
                    });

                    QRDecomposition decomp = new QRDecomposition(matrix);

                    DenseVector output = new DenseVector(expected.Length);

                    decomp.BasicSolution(inputOutput, ref output);

                    CheckEqual(expected, output);
                }

                public void x22()
                {
                    var matrix = new SquareMatrix(new Scalar[,] {
                        { 2, 1, },
                        {-1, 0, },
                    });

                    var inputOutput = new DenseVector(new Scalar[] {
                        3,
                        -4,
                    });

                    var expected = new DenseVector(new Scalar[] {
                         4,
                        -5,
                    });

                    QRDecomposition decomp = new QRDecomposition(matrix);

                    DenseVector output = new DenseVector(expected.Length);

                    decomp.BasicSolution(inputOutput, ref output);

                    CheckEqual(expected, output);
                }

                public void x33()
                {
                    var matrix = new SquareMatrix(new Scalar[,] {
                        { 2,  1,  3, },
                        {-1,  0,  7, },
                        { 0, -1, -1, },
                    });

                    var inputOutput = new DenseVector(new Scalar[] {
                        1,
                        2,
                        3,
                    });

                    var expected = new DenseVector(new Scalar[] {
                        1.5,
                        -3.5,
                        .5
                    });

                    QRDecomposition decomp = new QRDecomposition(matrix);

                    DenseVector output = new DenseVector(expected.Length);

                    decomp.BasicSolution(inputOutput, ref output);

                    CheckEqual(expected, output);
                }
            }

            public class RankDeficientSquare : TestFixture
            {
                public void x22R1()
                {
                    var matrix = new RectangularMatrix(new Scalar[,] {
                        {1, 2, },
                        {2, 4, },
                    });

                    var inputOutput = new DenseVector(new Scalar[] {
                        1,
                        -1,
                    });

                    var expected = new DenseVector(new Scalar[] {
                         0,
                        -.1,
                    });

                    QRDecomposition decomp = new QRDecomposition(matrix);

                    DenseVector output = new DenseVector(expected.Length);

                    decomp.BasicSolution(inputOutput, ref output);

                    CheckEqual(expected, output);
                }

                public void x33R1()
                {
                    var matrix = new SquareMatrix(new Scalar[,] {
                        {1, 2, 3, },
                        {2, 4, 6, },
                        {3, 6, 9, },
                    });

                    var inputOutput = new DenseVector(new Scalar[] {
                        1,
                        2,
                        3,
                    });

                    var expected = new DenseVector(new Scalar[] {
                         0,
                         0,
                         1.0/3.0,
                    });

                    QRDecomposition decomp = new QRDecomposition(matrix);

                    DenseVector output = new DenseVector(expected.Length);

                    decomp.BasicSolution(inputOutput, ref output);

                    CheckEqual(expected, output);
                }

                public void x33R2()
                {
                    var matrix = new SquareMatrix(new Scalar[,] {
                        { 13, 18, 48,  },
                        { 18, 25, 66,  },
                        { 48, 66, 180, },
                    });

                    var inputOutput = new DenseVector(new Scalar[] {
                        13,
                        18,
                        48,
                    });

                    var expected = new DenseVector(new Scalar[] {
                         0,
                         .5,
                         1.0/12.0,
                    });

                    QRDecomposition decomp = new QRDecomposition(matrix);

                    DenseVector output = new DenseVector(expected.Length);

                    decomp.BasicSolution(inputOutput, ref output);

                    CheckEqual(expected, output);
                }
            }

            public class FullRankTall : TestFixture
            {
                public void x42()
                {
                    var matrix = new RectangularMatrix(new Scalar[,] {
                        { 1, 2, },
                        { 2, 4, },
                        { 3, 6, },
                        { 4, 8, },
                    });

                    var inputOutput = new DenseVector(new Scalar[] {
                        1,
                        2,
                        3,
                        4,
                    });

                    var expected = new DenseVector(new Scalar[] {
                         0,
                         .5,
                    });

                    QRDecomposition decomp = new QRDecomposition(matrix);

                    DenseVector output = new DenseVector(expected.Length);

                    decomp.BasicSolution(inputOutput, ref output);

                    CheckEqual(expected, output);
                }

                public void x63()
                {
                    var matrix = new RectangularMatrix(new Scalar[,] {
                        { 1, 4, 1, },
                        { 2, 3, 0, },
                        { 3, 2, 0, },
                        { 4, 1, 1, },
                    });

                    var inputOutput = new DenseVector(new Scalar[] {
                        1,
                        2,
                        3,
                        4,
                    });

                    var expected = new DenseVector(new Scalar[] {
                         1,
                         0,
                         0,
                    });

                    QRDecomposition decomp = new QRDecomposition(matrix);

                    DenseVector output = new DenseVector(expected.Length);

                    decomp.BasicSolution(inputOutput, ref output);

                    CheckEqual(expected, output);
                }
            }

            public class RankDeficientTall : TestFixture
            {
                public void x42()
                {
                    var matrix = new RectangularMatrix(new Scalar[,] {
                        { 1, 1, },
                        { 2, 2, },
                        { 3, 3, },
                        { 4, 4, },
                    });

                    var inputOutput = new DenseVector(new Scalar[] {
                        1,
                        2,
                        3,
                        4,
                    });

                    var expected = new DenseVector(new Scalar[] {
                         1,
                         0,
                    });

                    QRDecomposition decomp = new QRDecomposition(matrix);

                    DenseVector output = new DenseVector(expected.Length);

                    decomp.BasicSolution(inputOutput, ref output);

                    CheckEqual(expected, output);
                }

                public void x63()
                {
                    var matrix = new RectangularMatrix(new Scalar[,] {
                        { 1, 2, 3, },
                        { 2, 4, 6, },
                        { 3, 6, 9, },
                        { 4, 8, 12, },
                    });

                    var inputOutput = new DenseVector(new Scalar[] {
                        1,
                        2,
                        3,
                        4,
                    });

                    var expected = new DenseVector(new Scalar[] {
                         0,
                         0,
                         1.0/3.0,
                    });

                    QRDecomposition decomp = new QRDecomposition(matrix);

                    DenseVector output = new DenseVector(expected.Length);

                    decomp.BasicSolution(inputOutput, ref output);

                    CheckEqual(expected, output);
                }
            }

            public class RankDeficientFat : TestFixture
            {
                public void x24()
                {
                    var matrix = new RectangularMatrix(new Scalar[,] {
                        { 1, 2, 3, 4, },
                        { 2, 4, 6, 8, },
                    });

                    var inputOutput = new DenseVector(new Scalar[] {
                        1,
                        2,
                    });

                    var expected = new DenseVector(new Scalar[] {
                         0,
                         0,
                         0,
                         .25,
                    });

                    QRDecomposition decomp = new QRDecomposition(matrix);

                    DenseVector output = new DenseVector(expected.Length);

                    decomp.BasicSolution(inputOutput, ref output);

                    CheckEqual(expected, output);
                }

                public void x36()
                {
                    var matrix = new RectangularMatrix(new Scalar[,] {
                        { 1, 2, 3, 4, },
                        { 2, 4, 6, 8, },
                        { 3, 6, 9, 12, },
                    });

                    var inputOutput = new DenseVector(new Scalar[] {
                        1,
                        2,
                        3,
                    });

                    var expected = new DenseVector(new Scalar[] {
                         0,
                         0,
                         0,
                         0.25,
                    });

                    QRDecomposition decomp = new QRDecomposition(matrix);

                    DenseVector output = new DenseVector(expected.Length);

                    decomp.BasicSolution(inputOutput, ref output);

                    CheckEqual(expected, output);
                }
            }

            public class FullRankFat : TestFixture
            {
                public void x24()
                {
                    var matrix = new RectangularMatrix(new Scalar[,] {
                        { 1, 2, 3, 4, },
                        { 1, 1, 1, 1, },
                    });

                    var inputOutput = new DenseVector(new Scalar[] {
                        1,
                        2,
                    });

                    var expected = new DenseVector(new Scalar[] {
                         7.0/3.0,
                         0,
                         0,
                         -1.0/3.0,
                    });

                    QRDecomposition decomp = new QRDecomposition(matrix);

                    DenseVector output = new DenseVector(expected.Length);

                    decomp.BasicSolution(inputOutput, ref output);

                    CheckEqual(expected, output);
                }

                public void x36()
                {
                    var matrix = new RectangularMatrix(new Scalar[,] {
                        { 1, 2, 3, 4, },
                        { 4, 3, 2, 1, },
                        { 1, 0, 0, 1, },
                    });

                    var inputOutput = new DenseVector(new Scalar[] {
                        1,
                        1,
                        1,
                    });

                    var expected = new DenseVector(new Scalar[] {
                         .6,
                         -.6,
                         0,
                         .4,
                    });

                    QRDecomposition decomp = new QRDecomposition(matrix);

                    DenseVector output = new DenseVector(expected.Length);

                    decomp.BasicSolution(inputOutput, ref output);

                    CheckEqual(expected, output);
                }
            }
        }

        public class ProjectToLeftNullSpace : TestFixture
        {
            public void IdentityMatrix()
            {
                var matrix = new RectangularMatrix(new Scalar[,]
                {
                    { 1, 0 },
                    {0, 1, },
                });

                var qr = new QRDecomposition(matrix);

                var expected = new DenseVector(new Scalar[] { 0, 0, });
                var vec = new DenseVector(new Scalar[] { 100, -100 });

                qr.ProjectToLeftNullSpace(ref vec);

                CheckEqual(expected, vec);
            }

            public void ZeroMatrix()
            {
                var matrix = new RectangularMatrix(new Scalar[,]
                {
                    { 0, 0 },
                    {0, 0, },
                });

                var qr = new QRDecomposition(matrix);

                var expected = new DenseVector(new Scalar[] { 100, -100, });
                var vec = new DenseVector(new Scalar[] { 100, -100 });

                qr.ProjectToLeftNullSpace(ref vec);

                CheckEqual(expected, vec);
            }

            public void RankDegenerate_Rows()
            {
                var matrix = new RectangularMatrix(new Scalar[,]
                {
                    { 1, 2 },
                    {3, 6, },
                });

                var qr = new QRDecomposition(matrix);

                var expected = new DenseVector(new Scalar[] { 120, -40, });
                var vec = new DenseVector(new Scalar[] { 100, -100 });

                qr.ProjectToLeftNullSpace(ref vec);

                CheckEqual(expected, vec);
            }

            public void RankDegenerate_Cols()
            {
                var matrix = new RectangularMatrix(new Scalar[,]
                {
                    { 1, 3, },
                    { 2, 6, },
                });

                var qr = new QRDecomposition(matrix);

                var expected = new DenseVector(new Scalar[] { 120, -60, });
                var vec = new DenseVector(new Scalar[] { 100, -100 });

                qr.ProjectToLeftNullSpace(ref vec);

                CheckEqual(expected, vec);
            }
        }

        public class NullQ : TestFixture
        {
            public class FullRankSquare : TestFixture
            {
                RectangularMatrix matrix = new RectangularMatrix(new Scalar[,]
                {
                    { 1, 2, 1, },
                    { 2, 3, 4, },
                    { 1, 1, 1, },
                });

                QRDecomposition qr;

                public override void FixtureSetup()
                {
                    qr = new QRDecomposition(matrix);
                }

                public void PremultipyByNullQ()
                {
                    DenseVector x = new Scalar[] { };
                    var y = qr.PremultiplyByNullQ(x);

                    DenseVector expected = new Scalar[] { 0, 0, 0, };

                    CheckEqual(expected, y);
                }

                public void PostmultipyByNullQ()
                {
                    DenseVector x = new Scalar[] { 5, 11, -3 };
                    var y = qr.PostmultiplyByNullQ(x);

                    DenseVector expected = new Scalar[] { };

                    CheckEqual(expected, y);
                }
            }

            public class ZeroMatrixSquare : TestFixture
            {
                RectangularMatrix matrix = new RectangularMatrix(new Scalar[,]
                {
                    { 0, 0, 0, },
                    { 0, 0, 0, },
                    { 0, 0, 0, },
                });

                QRDecomposition qr;

                public override void FixtureSetup()
                {
                    qr = new QRDecomposition(matrix);
                }

                public void PremultipyByNullQ()
                {
                    DenseVector x = new Scalar[] { 1, 2, 3, };
                    var y = qr.PremultiplyByNullQ(x);

                    DenseVector expected = new Scalar[] { 1, 2, 3, };

                    CheckEqual(expected, y);
                }

                public void PostmultipyByNullQ()
                {
                    DenseVector x = new Scalar[] { 5, 11, -3 };
                    var y = qr.PostmultiplyByNullQ(x);

                    DenseVector expected = new Scalar[] { 5, 11, -3 };

                    CheckEqual(expected, y);
                }
            }

            public class RankDeficientSquare : TestFixture
            {
                RectangularMatrix matrix = new RectangularMatrix(new Scalar[,]
                {
                    { 1, 2, 1, },
                    { 2, 3, 2, },
                    { 1, 1, 1, },
                });

                QRDecomposition qr;

                public override void FixtureSetup()
                {
                    qr = new QRDecomposition(matrix);
                }

                public void PremultipyByNullQ()
                {
                    DenseVector x = new Scalar[] { 7, };
                    var y = qr.PremultiplyByNullQ(x);

                    DenseVector expected = new Scalar[]
                    {
                        4.04145188432738,
                        -4.04145188432738,
                        4.04145188432738,
                    };

                    CheckEqual(expected, y);
                }

                public void PostmultipyByNullQ()
                {
                    DenseVector x = new Scalar[] { 5, 11, -3 };
                    var y = qr.PostmultiplyByNullQ(x);

                    DenseVector expected = new Scalar[] { -5.19615242270663 };

                    CheckEqual(expected, y);
                }
            }

            public class RankDeficientTall : TestFixture
            {
                RectangularMatrix matrix = new RectangularMatrix(new Scalar[,]
                {
                    { 3, 2, 1, },
                    { 5, 3, 2, },
                    { 2, 1, 1, },
                    { 5, 2, 3, },
                });

                QRDecomposition qr;

                public override void FixtureSetup()
                {
                    qr = new QRDecomposition(matrix);
                }

                public void PremultipyByNullQ()
                {
                    DenseVector x = new Scalar[] { 7, 11, };
                    var y = qr.PremultiplyByNullQ(x);

                    DenseVector expected = new Scalar[]
                    {
                        -9.86737753559965,
                        7.13974146633064,
                        3.77080281074534,
                        -2.72763606926898,
                    };

                    CheckEqual(expected, y);
                }

                public void PostmultipyByNullQ()
                {
                    DenseVector x = new Scalar[] { 5, 11, -3, 2, };
                    var y = qr.PostmultiplyByNullQ(x);

                    DenseVector expected = new Scalar[]
                    {
                        -4.31669152388046,
                        3.87722077709346,
                    };

                    CheckEqual(expected, y);
                }
            }

            public class FullRankFat : TestFixture
            {
                RectangularMatrix matrix = new RectangularMatrix(new Scalar[,]
                {
                    { 3, 2, 1, 0, },
                    { 5, 3, 2, 1, },
                    { 2, 1, 1, 0, },
                });

                QRDecomposition qr;

                public override void FixtureSetup()
                {
                    qr = new QRDecomposition(matrix);
                }

                public void PremultipyByNullQ()
                {
                    DenseVector x = new Scalar[] { };
                    var expectedX = x.Clone();

                    var y = qr.PremultiplyByNullQ(x);

                    DenseVector expected = new Scalar[]
                    {
                        0, 0, 0,
                    };

                    CheckEqual(expected, y);
                    CheckEqual(expectedX, x);
                }

                public void PostmultipyByNullQ()
                {
                    DenseVector x = new Scalar[] { 5, 11, -3, };
                    var expectedX = x.Clone();
                    var y = qr.PostmultiplyByNullQ(x);

                    DenseVector expected = new Scalar[]
                    {
                        
                    };

                    CheckEqual(expected, y);
                    CheckEqual(expectedX, x);
                }
            }

            public class RankDeficientFat : TestFixture
            {
                RectangularMatrix matrix = new RectangularMatrix(new Scalar[,]
                {
                    { 3, 2, 1, 2, },
                    { 5, 3, 2, 4, },
                    { 2, 1, 1, 2, },
                });

                QRDecomposition qr;

                public override void FixtureSetup()
                {
                    qr = new QRDecomposition(matrix);
                }

                public void PremultipyByNullQ()
                {
                    DenseVector x = new Scalar[] { 7, };
                    var y = qr.PremultiplyByNullQ(x);

                    DenseVector expected = new Scalar[]
                    {
                        4.04145188432738,
                        -4.04145188432738,
                        4.04145188432738,
                    };

                    CheckEqual(expected, y);
                }

                public void PostmultipyByNullQ()
                {
                    DenseVector x = new Scalar[] { 5, 11, -3, };
                    var y = qr.PostmultiplyByNullQ(x);

                    DenseVector expected = new Scalar[]
                    {
                        -5.19615242270663,
                    };

                    CheckEqual(expected, y);
                }
            }
        }

        public class AddColumnTests : TestFixture
        {
            public void EmptyCtor_Test()
            {
                var qr = new QRDecomposition(6, 6);

                CheckEqual(new RectangularMatrix(6, 0), qr.Q);
                CheckEqual(new RectangularMatrix(0, 0), qr.R);
                CheckEqual(0, qr.Rank);
                CheckEqual(new int[] { }, qr.P);
                CheckEqual(6, qr.RowCount);
                CheckEqual(0, qr.ColumnCount);
            }

            public void FromEmptyCtor_x11()
            {
                var qr = new QRDecomposition(1, 1);

                Check(qr.AddColumn(new Scalar[] { 7, }));

                CheckEqual(new RectangularMatrix(new Scalar[,] { { 1 } }), qr.Q);
                CheckEqual(new RectangularMatrix(new Scalar[,] { { 7 } }), qr.R);
                CheckEqual(1, qr.Rank);
                CheckEqual(new int[] { 0 }, qr.P);
                CheckEqual(1, qr.RowCount);
                CheckEqual(1, qr.ColumnCount);
            }

            public void FromEmptyCtor_x63()
            {
                var qr = new QRDecomposition(6, 3);

                Check(qr.AddColumn(new Scalar[] { 1, 2, 3, 4, 5, 6, }));
                Check(qr.AddColumn(new Scalar[] { 2, 2, 2, 2, 2, 2, }));
                Check(qr.AddColumn(new Scalar[] { 5, 2, 5, 2, 5, 2, }));

                RectangularMatrix expectedQ = new Scalar[,]
                {
                    { 0.1048284837,  0.7161148740,  0.2439750182, },
                    { 0.2096569673,  0.5012804118, -0.5367450401, },
                    { 0.3144854510,  0.2864459496,  0.3903600292, },
                    { 0.4193139347,  0.0716114874, -0.3903600292, },
                    { 0.5241424184, -0.1432229748,  0.5367450401, },
                    { 0.6289709020, -0.3580574370, -0.2439750182, },
                };

                RectangularMatrix expectedR = new Scalar[,]
                {
                    { 9.5393920142, 4.4027963142, 7.2331653734  },
                    { 0, 2.1483446221, 4.7263581687  },
                    { 0, 0, 3.5132402626  },
                };

                CheckEqual(expectedQ, qr.Q);
                CheckEqual(expectedR, qr.R);
                CheckEqual(3, qr.Rank);
                CheckEqual(new int[] { 0, 1, 2, }, qr.P);
                CheckEqual(6, qr.RowCount);
                CheckEqual(3, qr.ColumnCount);
            }

            public void FromInitialDecomposition_x63()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { 1, 4, },
                    { 2, 3, },
                    { 3, 2, },
                    { 4, 1, },
                });

                var expectedQ = new RectangularMatrix(new Scalar[,] {
                    { 1.0/Math.Sqrt(30),  2.0/Math.Sqrt(6),  0.5, },
                    { 2.0/Math.Sqrt(30),  1.0/Math.Sqrt(6), -0.5, },
                    { 3.0/Math.Sqrt(30),                 0, -0.5, },
                    { 4.0/Math.Sqrt(30), -1.0/Math.Sqrt(6),  0.5, },
                });

                var expectedR = new RectangularMatrix(new Scalar[,] {
                    { Math.Sqrt(30), 1.0 / Math.Sqrt(.075), 1.0/Math.Sqrt(1.2), },
                    {             0, 30.0 / Math.Sqrt(54),  1.0/Math.Sqrt(6), },
                    {             0,                   0,                  1, },
                });

                var expectedP = new int[] { 0, 1, 2, };

                QRDecomposition decomp = new QRDecomposition(matrix);
                Check(decomp.AddColumn(new Scalar[] { 1, 0, 0, 1, }));

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedP, decomp.P);
            }
        }

        public class RemoveColumnsFromEndTests : TestFixture
        {
            public void ZeroCount()
            {
                var matrix = new RectangularMatrix(new Scalar[,] 
                {
                   { 5, 1, 3, },
                   { 1, 5, 3, },
                   { 1, 0, 0, },
                });

                var qr = new QRDecomposition(matrix);

                var expectedQ = new RectangularMatrix(new Scalar[,]
                {
                    { 0.9622504486, -0.1804046420,  0.2037847865, },
                    { 0.1924500897,  0.9804600111, -0.0407569573, },
                    { 0.1924500897, -0.0784368009, -0.9781669751, },
                });

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { 5.1961524227, 1.9245008973, 3.4641016151, },
                    { 0, 4.7218954135, 2.4001661072, },
                    { 0, 0, 0.4890834876, },
                });

                var expectedP = new int[] 
                {
                    0, 1, 2,
                };

                qr.RemoveColumnsFromEnd(0);

                CheckEqual(expectedQ, qr.Q);
                CheckEqual(expectedR, qr.R);
                CheckEqual(expectedP, qr.P);
                CheckEqual(3, qr.RowCount);
                CheckEqual(3, qr.ColumnCount);
                CheckEqual(3, qr.Rank);
            }

            public void DefaultCount_1()
            {
                var matrix = new RectangularMatrix(new Scalar[,] 
                {
                   { 5, 1, 3, },
                   { 1, 5, 3, },
                   { 1, 0, 0, },
                });

                var qr = new QRDecomposition(matrix);

                var expectedQ = new RectangularMatrix(new Scalar[,]
                {
                    { 0.9622504486, -0.1804046420, },
                    { 0.1924500897,  0.9804600111, },
                    { 0.1924500897, -0.0784368009, },
                });

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { 5.1961524227, 1.9245008973, },
                    { 0, 4.7218954135, },
                });

                var expectedP = new int[] 
                {
                    0, 1,
                };

                qr.RemoveColumnsFromEnd();

                CheckEqual(expectedQ, qr.Q);
                CheckEqual(expectedR, qr.R);
                CheckEqual(expectedP, qr.P);
                CheckEqual(3, qr.RowCount);
                CheckEqual(2, qr.ColumnCount);
                CheckEqual(2, qr.Rank);
            }

            public void Multiple()
            {
                var matrix = new RectangularMatrix(new Scalar[,] 
                {
                   { 5, 1, 3, },
                   { 1, 5, 3, },
                   { 1, 0, 0, },
                });

                var qr = new QRDecomposition(matrix);

                var expectedQ = new RectangularMatrix(new Scalar[,]
                {
                    { 0.9622504486, },
                    { 0.1924500897, },
                    { 0.1924500897, },
                });

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { 5.1961524227, },
                });

                var expectedP = new int[] 
                {
                    0,
                };

                qr.RemoveColumnsFromEnd(2);

                CheckEqual(expectedQ, qr.Q);
                CheckEqual(expectedR, qr.R);
                CheckEqual(expectedP, qr.P);
                CheckEqual(3, qr.RowCount);
                CheckEqual(1, qr.ColumnCount);
                CheckEqual(1, qr.Rank);
            }

            public void All()
            {
                var matrix = new RectangularMatrix(new Scalar[,] 
                {
                   { 5, 1, 3, },
                   { 1, 5, 3, },
                   { 1, 0, 0, },
                });

                var qr = new QRDecomposition(matrix);

                var expectedQ = new RectangularMatrix(new Scalar[,]
                {
                    { },
                    { },
                    { },
                });

                var expectedR = new RectangularMatrix(0, 0);

                var expectedP = new int[] 
                {
                    
                };

                qr.RemoveColumnsFromEnd(3);

                CheckEqual(expectedQ, qr.Q);
                CheckEqual(expectedR, qr.R);
                CheckEqual(expectedP, qr.P);
                CheckEqual(3, qr.RowCount);
                CheckEqual(0, qr.ColumnCount);
                CheckEqual(0, qr.Rank);
            }

            public void TooMany()
            {
                var matrix = new RectangularMatrix(new Scalar[,] 
                {
                   { 5, 1, 3, },
                   { 1, 5, 3, },
                   { 1, 0, 0, },
                });

                var qr = new QRDecomposition(matrix);

                CheckThrow(typeof(Exception));
                qr.RemoveColumnsFromEnd(4);
            }

            public void Negative()
            {
                var matrix = new RectangularMatrix(new Scalar[,] 
                {
                   { 5, 1, 3, },
                   { 1, 5, 3, },
                   { 1, 0, 0, },
                });

                var qr = new QRDecomposition(matrix);

                CheckThrow(typeof(Exception));
                qr.RemoveColumnsFromEnd(-1);
            }
        }
    }
}