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 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 = inputVector.PostMatrixMultiplyBy(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(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.PremultiplyByQ(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 }); 

                outputVector1 = decomp.PremultiplyByQ(outputVector1);

                var outputVector2 = inputVector.PreMatrixMultiplyBy(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 expected = new DenseVector(new Scalar[] { 1, 2, 3, 4 });
                var inputVector = new DenseVector(new Scalar[] { Math.Sqrt(30.0) });
                var outputVector = inputVector.Clone();

                outputVector = decomp.PremultiplyByQ(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 = inputVector.Clone();

                outputVector = decomp.PremultiplyByQ(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 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 NotAnythingImportant : TestFixture
        {
            public void IdentityFirst()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    {  1, 0,  7, 9, },
                    {  0, 1, 10, 8, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);
            }

            public void IdentityLast()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    {  7, 9, 1, 0, },
                    {  10, 8, 0, 1, },
                });

                QRDecomposition decomp = new QRDecomposition(matrix);
            }
        }

        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);

                    decomp.BasicSolution(ref inputOutput);

                    CheckEqual(expected, inputOutput);
                }

                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);

                    decomp.BasicSolution(ref inputOutput);

                    CheckEqual(expected, inputOutput);
                }

                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);

                    decomp.BasicSolution(ref inputOutput);

                    CheckEqual(expected, inputOutput);
                }
            }

            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);

                    decomp.BasicSolution(ref inputOutput);

                    CheckEqual(expected, inputOutput);                    
                }

                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);

                    decomp.BasicSolution(ref inputOutput);

                    CheckEqual(expected, inputOutput);
                }

                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);

                    decomp.BasicSolution(ref inputOutput);

                    CheckEqual(expected, inputOutput);
                }
            }

            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);

                    decomp.BasicSolution(ref inputOutput);

                    CheckEqual(expected, inputOutput);
                }

                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);

                    decomp.BasicSolution(ref inputOutput);

                    CheckEqual(expected, inputOutput);
                }
            }

            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);

                    decomp.BasicSolution(ref inputOutput);

                    CheckEqual(expected, inputOutput);
                }

                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);

                    decomp.BasicSolution(ref inputOutput);

                    CheckEqual(expected, inputOutput);
                }
            }

            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);

                    decomp.BasicSolution(ref inputOutput);

                    CheckEqual(expected, inputOutput);
                }

                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);

                    decomp.BasicSolution(ref inputOutput);

                    CheckEqual(expected, inputOutput);
                }
            }

            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);

                    decomp.BasicSolution(ref inputOutput);

                    CheckEqual(expected, inputOutput);
                }

                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);

                    decomp.BasicSolution(ref inputOutput);

                    CheckEqual(expected, inputOutput);
                }
            }
        }
    }
}