using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using Azimuth.DenseLinearAlgebra;

namespace Azimuth.UnitTests.DenseLinearAlgebra
{
    public class CompleteOrthogonalDecompositionTests : UnitTestSharp.TestFixture
    {
        public class Ctor : UnitTestSharp.TestFixture
        {
            public void IdentityMatrix()
            {
                var matrix = new RectangularMatrix(new Scalar[,]
                {
                    { 1, 0, },
                    { 0, 1, },
                });

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

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

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

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(4, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }

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

                var expectedQ = new RectangularMatrix(2, 0);

                var expectedR = new RectangularMatrix(0, 0);

                var expectedK = new RectangularMatrix(0, 2);

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(4, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }

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

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

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

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

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(4, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }

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

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

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { 3, 2, },
                    { 0, 1, },
                });

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

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(4, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }

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

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

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { 5, 2.4, },
                    { 0, 1.8, },
                });

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

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(4, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }

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

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

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { 5, 2.4, },
                    { 0, 1.8, },
                });

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

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(4, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }

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

                var expectedQ = new RectangularMatrix(new Scalar[,] 
                {
                    { 0.4242640687, 0.3352576037, -0.8411910242, },
                    { 0.5656854249, 0.6272561617, 0.535303379, },
                    { 0.7071067812, -0.7029594916, 0.0764719113, },
                });

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { 7.0710678119, 3.2526911935, 2.2627416998, },
                    { 0, 1.8493242009, 0.8868104355, },
                    { 0, 0, 0.3058876452, },
                });

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

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(6, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }

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

                var expectedQ = new RectangularMatrix(new Scalar[,] 
                {
                    { 0.3487429162, -0.7382716607 },
                    { 0.464990555, 0.6711560552 },
                    { 0.8137334712, -0.0671156055 },
                });

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { -9.4233751915,  -6.0825403654 }, 
                    { 0,  -0.4502251689 },
                });

                var expectedK = new RectangularMatrix(new Scalar[,]
                {
                    { 0.1825741858, -0.3651483717, -0.9128709292 },
                    { -0.894427191, -0.4472135955, 0 },
                });

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(6, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }

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

                var expectedQ = new RectangularMatrix(new Scalar[,] 
                {
                    { 0.4866642634, 0.6556100681 },
                    { 0.8111071057, -0.0936585812 },
                    { 0.3244428423, -0.7492686493 },
                });

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { -7.5498344353, -0.9176629355 }, 
                    { 0, -0.3973597071 },
                });

                var expectedK = new RectangularMatrix(new Scalar[,]
                {
                    { -0.4082482905, -0.4082482905, -0.8164965809, },
                    { 0.7071067812, -0.7071067812, 0, }, 
                });

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(6, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }

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

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

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { 6.164414003, 0, 0 },
                    { 3.731092686, 4.0098562778, 0 },
                    { 2.433321317, 1.9753956458, 0.4204276523 },
                });

                var expectedK = new RectangularMatrix(new Scalar[,]
                {
                    { 0.3244428423, 0.4866642634, 0.8111071057, 0},
                    { -0.0525022099, 0.0459394336, -0.0065627762, 0.9975419873},
                    { 0.7474269374, -0.6539985702, 0.0934283672, 0.0700712754},
                });

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(matrix.ColumnCount + matrix.RowCount, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }

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

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

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { 5.4772255751, 0, 0 },
                    { 4.7469288317, 2.3380903889, 0 },
                    { 3.8340579025, 0.3421595691, 0.4276994614 },
                });

                var expectedK = new RectangularMatrix(new Scalar[,]
                {
                    { 0.3651483717, 0.7302967433, 0.5477225575, 0.1825741858 },
                    { -0.3136462717, 0.6558058408, -0.6843191382, 0.0570265949 },
                    { -0.6843191382, -0.0570265949, 0.3136462717, 0.6558058408 },
                });

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                decomp.R.ToString();

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(matrix.ColumnCount + matrix.RowCount, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }

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

                var expectedQ = new RectangularMatrix(new Scalar[,] 
                {
                    { -0.4082482905, 0.7071067812 },
                    { -0.4082482905, -0.7071067812 },
                    { -0.8164965809, 0 },
                });

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { -13.0766968306, 0 },
                    { -0.5298129428, -3.2740339408 },
                });

                var expectedK = new RectangularMatrix(new Scalar[,]
                {
                    { 0.2809757435, 0.4682929058, 0.7492686493, 0.3746343246 },
                    { 0.1705059281, 0.1401937631, 0.3106996912, -0.9245210323 },
                });

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                decomp.R.ToString();

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(matrix.ColumnCount + matrix.RowCount, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }

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

                var expectedQ = new RectangularMatrix(new Scalar[,] 
                {
                    { 0.3015113446, 0.4487745552, -0.4837438342 },
                    { 0.4020151261, 0.7354916321, 0.388892102 },
                    { 0.5025189076, -0.2119213177, -0.644991779 },
                    { 0.7035264707, -0.4612405151, 0.4458031414 },
                });

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { 9.9498743711, 3.0151134458, 3.0151134458 },
                    { 0, 2.430862174, 0.7853554716 },
                    { 0, 0, 0.5406548736 },
                });

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

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(7, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }

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

                var expectedQ = new RectangularMatrix(new Scalar[,] 
                {
                    { 0.4200840252, 0.3055661657, -0.6299407883 },
                    { 0.5601120336, 0.5805757148, 0.1259881577 },
                    { 0.700140042, -0.7028021811, 0.1259881577 },
                    { 0.1400280084, 0.2750095491, 0.755928946 },
                });

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { 7.1414284285, 3.3606722017, 2.3804761428 },
                    { 0, 1.9250668438, 1.0389249633 },
                    { 0, 0, 0.5039526307 },
                });

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

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(7, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }

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

                var expectedQ = new RectangularMatrix(new Scalar[,] 
                {
                    { 0.2153874476, 0.2846324908, -0.7890453679 },
                    {  0.2871832634, 0.5582112927, 0.5965952782 },
                    {  0.3589790793, -0.7765216497, 0.1347150628 },
                    {  0.8615497903, 0.0663221338, -0.0577350269 },
                });

                var expectedR = new RectangularMatrix(new Scalar[,]
                {
                    { 13.928388277, 6.8206025069, 4.5949322152 },
                    { 0, 1.8653100127, 0.8898219616 },
                    { 0, 0, 0.3079201436 },
                });

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

                int progress = 0;

                var decomp = new CompleteOrthogonalDecomposition(ref matrix, ref progress);

                CheckEqual(expectedQ, decomp.Q);
                CheckEqual(expectedR, decomp.R);
                CheckEqual(expectedK, decomp.K);
                CheckEqual(7, progress);
                CheckEqual(matrix.RowCount, decomp.RowCount);
                CheckEqual(matrix.ColumnCount, decomp.ColumnCount);
            }
        }

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

                CompleteOrthogonalDecomposition decomp;

                public override void FixtureSetup()
                {
                    decomp = new CompleteOrthogonalDecomposition(ref matrix);
                }

                public void PremultiplyByQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        9.5372810542,
                        9.0074226688,
                        2.6280424181,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PostmultiplyByQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        6.5053823864,
                        10.1875052474,
                        -5.7353933471,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        93.762359186,
                        10.2848380998,
                        -0.9176629355,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PostmultiplyByR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        77.7817459309000,
                        48.7248725348000,
                        30.1801688107000,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByInverseR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByInverseR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        0.7894752008,
                        8.488189036672765,
                        -9.807522623,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByInverseR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PostmultiplyByInverseR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        1.55563491860281,
                        1.04903185661845,
                        -24.3563037459,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        7,
                        -3,
                        11,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PostmultiplyByK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -3,
                        11,
                        7,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByNullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByNullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByNullK()
                {
                    var vector = new DenseVector(3);
                    vector.SetLength(0);
                    
                    decomp.PostmultiplyByNullK(ref vector);

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

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByNullQ()
                {
                    var vector = new DenseVector(3);
                    vector.SetLength(0);                    

                    decomp.PremultiplyByNullQ(ref vector);

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

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByNullQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, });
                    
                    decomp.PostmultiplyByNullQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByFullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });

                    decomp.PremultiplyByFullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        7,
                        -3,
                        11,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByFullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });

                    decomp.PostmultiplyByFullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -3,
                        11,
                        7,
                    };

                    CheckEqual(expectedOutput, vector);
                }
            }

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

                CompleteOrthogonalDecomposition decomp;

                public override void FixtureSetup()
                {
                    decomp = new CompleteOrthogonalDecomposition(ref matrix);
                }

                public void PremultiplyByQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        9.94257737406546,
                        8.26656809408475,
                        -1.67600927998071,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PostmultiplyByQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        10.0577281101073,
                        8.80390662886792,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7 });
                    
                    decomp.PremultiplyByR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -89.471819336354,
                        -2.78151794983659,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7 });
                    
                    decomp.PostmultiplyByR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -83.0481787879782,
                        -12.8758102401413,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByInverseR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7 });
                    
                    decomp.PremultiplyByInverseR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        0.684227923353702,
                        -17.6162803489651,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByInverseR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7 });
                    
                    decomp.PostmultiplyByInverseR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -1.45698559277155,
                        -14.2515162521968,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -8.98146239020499,
                        9.89949493661166,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -1 });
                    vector.SetLength(2);

                    decomp.PostmultiplyByK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        0.459016273203332,
                        -8.98146239020499,
                        -9.44047866340832,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByNullK()
                {
                    var vector = new DenseVector(new Scalar[] { -1, -1, -1, });
                    vector.SetLength(1);
                    vector[0] = 7;
                    
                    decomp.PostmultiplyByNullK(ref vector);

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

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByNullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, });
                    
                    decomp.PremultiplyByNullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        0.577350269189625,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByNullQ()
                {
                    var vector = new DenseVector(new Scalar[] { 7, -1, -1, });
                    vector.SetLength(1);

                    decomp.PremultiplyByNullQ(ref vector);

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

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByNullQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, });
                    
                    decomp.PostmultiplyByNullQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        0.577350269189625,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByFullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });

                    decomp.PremultiplyByFullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -8.98146239020499,
                        9.89949493661166,
                        0.57735026918966,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByFullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -1 });
                    
                    decomp.PostmultiplyByFullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -0.118333996,
                        -8.404112121,
                        -10.0178289326
                    };

                    CheckEqual(expectedOutput, vector);
                }
            }

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

                CompleteOrthogonalDecomposition decomp;

                public override void FixtureSetup()
                {
                    decomp = new CompleteOrthogonalDecomposition(ref matrix);
                }

                public void PremultiplyByQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, -1 });
                    vector.SetLength(3);

                    decomp.PremultiplyByQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        8.64970980207106,
                        9.84729790020575,
                        2.40396072176009,
                        1.19758809813469,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, 2, });
                    decomp.PostmultiplyByQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        6.72134440333445,
                        10.0836834673103,
                        -4.91353814911996,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        94.9389896970991,
                        10.3586930164188,
                        -1.51185789203691,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PostmultiplyByR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        78.5557127139714,
                        50.4428621246542,
                        31.9458544223774,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByInverseR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByInverseR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        0.30159315590635,
                        6.84893539212367,
                        -5.95294044989533,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByInverseR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PostmultiplyByInverseR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        1.54030809243081,
                        0.947255113595818,
                        -15.181572999204,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        7,
                        -3,
                        11,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PostmultiplyByK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -3,
                        11,
                        7,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByNullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByNullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByNullK()
                {
                    var vector = new DenseVector(new Scalar[] { -1, -1, -1});
                    vector.SetLength(0);

                    decomp.PostmultiplyByNullK(ref vector);

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

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByNullQ()
                {
                    var vector = new DenseVector(new Scalar[] { 7, -1, -1, -1, });
                    vector.SetLength(1);

                    decomp.PremultiplyByNullQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        4.041451884327380,
                        -4.041451884327380,
                        0,
                        4.041451884327380,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByNullQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, 2, });
                    
                    decomp.PostmultiplyByNullQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        3.46410161513775,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByFullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });

                    decomp.PremultiplyByFullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        7,
                        -3,
                        11,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByFullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });

                    decomp.PostmultiplyByFullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -3,
                        11,
                        7,
                    };

                    CheckEqual(expectedOutput, vector);
                }
            }

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

                CompleteOrthogonalDecomposition decomp;

                public override void FixtureSetup()
                {
                    decomp = new CompleteOrthogonalDecomposition(ref matrix);
                }

                public void PremultiplyByQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -1, -1, });
                    vector.SetLength(2);

                    decomp.PremultiplyByQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        8.00995227739239,
                        6.19605592783456,
                        -1.81389634955782,
                        8.00995227739239,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, 2, });
                    
                    decomp.PostmultiplyByQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        9.91881942185683,
                        6.26234950131305,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7 });
                    
                    decomp.PremultiplyByR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -100.302649526799,
                        -3.2288592281011,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7 });
                    
                    decomp.PostmultiplyByR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -92.3607059306067,
                        -15.7090563078325,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByInverseR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7 });
                    
                    decomp.PremultiplyByInverseR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        0.740520083785624,
                        -15.1756383720752,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByInverseR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7 });
                    
                    decomp.PostmultiplyByInverseR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -1.31008093518591,
                        -11.9532653422627,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -8.98146239020499,
                        9.89949493661166,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, 0 });
                    
                    vector.SetLength(2);

                    decomp.PostmultiplyByK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        0.459016273203338,
                        -8.98146239020499,
                        -9.44047866340832,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByNullQ()
                {
                    var vector = new DenseVector(new Scalar[] { 7, 11, -1, -1, });
                    vector.SetLength(2);

                    decomp.PremultiplyByNullQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -5.12542376137919,
                        -4.77783712229788,
                        4.77783712183038,
                        9.90326088398874
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByNullQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, 2, });
                    
                    decomp.PostmultiplyByNullQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -1.8846098017,
                        -6.46902202000333,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByNullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByNullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        0.5773502692,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByNullK()
                {
                    var vector = new DenseVector(new Scalar[] { 7, -1, -1, });
                    vector.SetLength(1);
                    
                    decomp.PostmultiplyByNullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        4.04145188416239,
                        -4.04145188416239,
                        4.04145188416239,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByFullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });

                    decomp.PremultiplyByFullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -8.98146239020499,
                        9.89949493661166,
                        0.57735026918966,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByFullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -1 });

                    decomp.PostmultiplyByFullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -0.118333996,
                        -8.404112121,
                        -10.0178289326
                    };

                    CheckEqual(expectedOutput, vector);
                }
            }

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

                CompleteOrthogonalDecomposition decomp;

                public override void FixtureSetup()
                {
                    decomp = new CompleteOrthogonalDecomposition(ref matrix);
                }

                public void PremultiplyByQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -3,
                        11, 
                        7,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PostmultiplyByQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        7,
                        -3,
                        11,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        60.2494813255683,
                        68.582849871126,
                        43.2866555274966,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PostmultiplyByR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        81.9758094399399,
                        15.3401540149782,
                        -1.28309838415245,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByInverseR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PremultiplyByInverseR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        2.00831604418561,
                        -1.08350530217318,
                        -24.1507629194917,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByInverseR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PostmultiplyByInverseR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        3.43398091546612,
                        4.02037493701102,
                        -7.01427116670008,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, 2 });
                    
                    decomp.PremultiplyByK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        7.85068999090738,
                        3.30754250137077,
                        -7.55602381778666,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, 0 });
                    vector.SetLength(3);

                    decomp.PostmultiplyByK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        3.87406560124317,
                        12.7949848468193,
                        0.293775350009208,
                        0.440084685777055,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByNullQ()
                {
                    var vector = new DenseVector(new Scalar[] { -1, -1, -1, });
                    vector.SetLength(0);

                    decomp.PremultiplyByNullQ(ref vector);

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

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByNullQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, });
                    
                    decomp.PostmultiplyByNullQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByNullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, 2 });
                    
                    decomp.PremultiplyByNullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        7.30296743347978,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByNullK()
                {
                    var vector = new DenseVector(new Scalar[] { 7, -1, -1, -1, });
                    vector.SetLength(1);

                    decomp.PostmultiplyByNullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        3.83405790266370,
                        -1.27801930092191,
                        -2.55603860164340,
                        5.1120772034,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByFullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, 2 });

                    decomp.PremultiplyByFullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        7.85068999090738,
                        3.30754250137077,
                        -7.55602381778666,
                        7.3029674334,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByFullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, 2 });
                    
                    decomp.PostmultiplyByFullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        4.9695107163,
                        12.429836475,
                        -0.4365213933,
                        1.9006781725,
                    };

                    CheckEqual(expectedOutput, vector);
                }
            }

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

                CompleteOrthogonalDecomposition decomp;

                public override void FixtureSetup()
                {
                    decomp = new CompleteOrthogonalDecomposition(ref matrix);
                }

                public void PremultiplyByQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -1 });
                    vector.SetLength(2);

                    decomp.PremultiplyByQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -9.44047866340832,
                        -8.98146239020499,
                        0.459016273203338,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3 });
                    
                    decomp.PostmultiplyByQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -8.98146239020499,
                        -9.89949493661167,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7 });
                    
                    decomp.PremultiplyByR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -121.249742267767,
                        2.78631180791942,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7 });
                    
                    decomp.PostmultiplyByR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -114.10010703577,
                        -8.4488292709335,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByInverseR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7 });
                    
                    decomp.PremultiplyByInverseR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -0.997940265578332,
                        -6.64410500935656,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByInverseR()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7 });
                    
                    decomp.PostmultiplyByInverseR(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -1.53533975795832,
                        -5.79962009275944,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, 2 });
                    
                    decomp.PremultiplyByK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        8.11111111111111,
                        -3.09560299443292,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, 0, 0 });
                    vector.SetLength(2);

                    decomp.PostmultiplyByK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        3.76079898880611,
                        6.84129456343253,
                        9.74214917910511,
                        3.76079898880611,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByNullQ()
                {
                    var vector = new DenseVector(new Scalar[] { 7, -1, -1, });
                    vector.SetLength(1);

                    decomp.PremultiplyByNullQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        4.04145188416239,
                        -4.04145188416239,
                        4.04145188416239,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByNullQ()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, });
                    
                    decomp.PostmultiplyByNullQ(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        0.5773502692,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByNullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, 2 });
                    
                    decomp.PremultiplyByNullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -10.20792698365459,
                        -1.85076885182660,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByNullK()
                {
                    var vector = new DenseVector(new Scalar[] { 7, 11, -1, -1});
                    vector.SetLength(2);
                    
                    decomp.PostmultiplyByNullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        -9.904789748728298,
                        0.205449116673616,
                        0.410898233376908,
                        8.466645931876013,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PremultiplyByFullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, 2 });

                    decomp.PremultiplyByFullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        8.11111111111111,
                        -3.09560299443292,
                        0,
                        -1.850768852,
                    };

                    CheckEqual(expectedOutput, vector);
                }

                public void PostmultiplyByFullK()
                {
                    var vector = new DenseVector(new Scalar[] { 11, 7, -3, -2 });
                    
                    decomp.PostmultiplyByFullK(ref vector);

                    DenseVector expectedOutput = new Scalar[] 
                    {
                        4.4181864225,
                        7.011675271,
                        10.082910594,
                        1.9107466019,
                    };

                    CheckEqual(expectedOutput, vector);
                }
            }
        }

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

                var qr = new CompleteOrthogonalDecomposition(ref 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 CompleteOrthogonalDecomposition(ref 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 CompleteOrthogonalDecomposition(ref 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 CompleteOrthogonalDecomposition(ref 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 ProjectToColumnSpace : UnitTestSharp.TestFixture
        {
            public void Square_FullRank()
            {
                var matrix = new RectangularMatrix(new Scalar[,] 
                {
                    { 1, 2, 3, },
                    { 3, 2, 1, },
                    { 1, 4, 4, },
                });

                var cod = new CompleteOrthogonalDecomposition(ref matrix);

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

                var actual = input.Clone();

                cod.ProjectToColumnSpace(ref actual);

                CheckEqual(expected, actual);
            }

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

                var cod = new CompleteOrthogonalDecomposition(ref matrix);

                DenseVector input = new Scalar[] { 1, 5, 7, 2 };
                DenseVector expected = new Scalar[] { 2.94871794871795, 5.27838827838828, 5.88644688644689, 1.58241758241758};

                var actual = input.Clone();
                cod.ProjectToColumnSpace(ref actual);

                CheckEqual(expected, actual);
            }

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

                var cod = new CompleteOrthogonalDecomposition(ref matrix);

                DenseVector input = new Scalar[] { 1, 5, 7, };
                DenseVector expected = new Scalar[] { 1, 5, 7, };

                var actual = input.Clone();
                cod.ProjectToColumnSpace(ref actual);

                CheckEqual(expected, actual);
            }

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

                var cod = new CompleteOrthogonalDecomposition(ref matrix);

                DenseVector input = new Scalar[] { 1, 5, 13, };
                DenseVector expected = new Scalar[] { 3.55555555555556, 7.55555555555556, 2.77777777777778 };

                var actual = input.Clone();
                cod.ProjectToColumnSpace(ref actual);

                CheckEqual(expected, actual);
            }
        }

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

                var qr = new CompleteOrthogonalDecomposition(ref matrix);

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

                qr.ProjectToNullSpace(ref vec);

                CheckEqual(expected, vec);
            }

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

                var qr = new CompleteOrthogonalDecomposition(ref matrix);

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

                qr.ProjectToNullSpace(ref vec);

                CheckEqual(expected, vec);
            }

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

                var qr = new CompleteOrthogonalDecomposition(ref matrix);

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

                var shouldBeZero = new DenseVector(2);
                shouldBeZero.Assignment(vec);

                qr.PremultiplyByMatrix(ref shouldBeZero);

                CheckEqual(expected, vec);
                CheckEqual(new DenseVector(2), shouldBeZero);
            }

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

                var qr = new CompleteOrthogonalDecomposition(ref matrix);

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

                var shouldBeZero = vec.Clone();
                qr.PremultiplyByMatrix(ref shouldBeZero);

                CheckEqual(expected, vec);
                CheckEqual(new DenseVector(2), shouldBeZero);
            }

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

                var qr = new CompleteOrthogonalDecomposition(ref matrix);

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

                var shouldBeZero = vec.Clone();
                qr.PremultiplyByMatrix(ref shouldBeZero);

                CheckEqual(expected, vec);
                CheckEqual(new DenseVector(2), shouldBeZero);
            }

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

                var qr = new CompleteOrthogonalDecomposition(ref matrix);

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

                var shouldBeZero = vec.Clone();
                qr.PremultiplyByMatrix(ref shouldBeZero);

                CheckEqual(expected, vec);
                CheckEqual(new DenseVector(2), shouldBeZero);
            }

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

                var qr = new CompleteOrthogonalDecomposition(ref matrix);

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

                var shouldBeZero = new DenseVector(3);
                shouldBeZero.SetLength(2);
                shouldBeZero.Assignment(vec);
                qr.PremultiplyByMatrix(ref shouldBeZero);

                CheckEqual(expected, vec);
                CheckEqual(new DenseVector(3), shouldBeZero);
            }

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

                var qr = new CompleteOrthogonalDecomposition(ref matrix);

                var vec = new DenseVector(new Scalar[] { 100, -100, 100, });
                CheckThrow(typeof(Exception));
                qr.ProjectToNullSpace(ref vec);
            }

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

                var qr = new CompleteOrthogonalDecomposition(ref matrix);

                var expected = new DenseVector(new Scalar[] 
                {
                    19.73684210526322,
                    7.89473684210527,
                    -11.84210526315793,
                });
                var vec = new DenseVector(new Scalar[] { 100, -100, 50, });
                qr.ProjectToNullSpace(ref vec);

                var shouldBeZero = vec.Clone();
                qr.PremultiplyByMatrix(ref shouldBeZero);

                CheckEqual(expected, vec);
                CheckEqual(new DenseVector(2), shouldBeZero);
            }

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

                var qr = new CompleteOrthogonalDecomposition(ref matrix);

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

                CheckThrow(typeof(Exception));
                qr.ProjectToNullSpace(ref vec);
            }
        }

        public class UnconstrainedLeastSquaresSolution : UnitTestSharp.TestFixture
        {
            public class FullRankSquare : UnitTestSharp.TestFixture
            {
                public void x11()
                {
                    var matrix = new RectangularMatrix(new Scalar[,] {
                        { 7 },
                    });

                    var input = new DenseVector(new Scalar[] {
                        28,
                    });

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

                    
                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    var output = input.Clone();
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }

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

                    var input = new DenseVector(new Scalar[] {
                        4,
                        2,
                    });

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

                    DenseVector output = input.Clone();

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }

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

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

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

                    DenseVector output = input.Clone();

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }

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

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

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

                    DenseVector output = input.Clone();

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }

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

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

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

                    DenseVector output = input.Clone();

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }

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

                    var input = new DenseVector(new Scalar[] {
                         4,
                         2,
                         0,
                    });

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

                    DenseVector output = input.Clone();

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }
            }

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

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

                    var expected = new DenseVector(new Scalar[] {
                        .12,
                        .24,
                    });

                    DenseVector output = input.Clone();

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }

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

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

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

                    DenseVector output = input.Clone();

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }
            }

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

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

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

                    DenseVector output = new DenseVector(4);
                    output.SetLength(2);
                    output.Assignment(input);

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }

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

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

                    var expected = new DenseVector(new Scalar[] {
                        -6.0/28.0,
                         5.0/28.0,
                         5.0/28.0,
                         5.0/28.0,
                         5.0/28.0,
                        -6.0/28.0,
                    });

                    DenseVector output = new DenseVector(expected.Length);
                    output.SetLength(input.Length);
                    output.Assignment(input);

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }
            }

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

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

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

                    DenseVector output = input.Clone();

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }

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

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

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

                    DenseVector output = input.Clone();

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }
            }

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

                    var input = new DenseVector(new Scalar[] {
                        100,
                        100,
                    });

                    var expected = new DenseVector(new Scalar[] {
                        2,
                        4,
                        6,
                        8,
                    });

                    DenseVector output = new DenseVector(4);
                    output.SetLength(2);
                    output.Assignment(input);

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }

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

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

                    var expected = new DenseVector(new Scalar[] {
                        1 * 3.0/637.0,
                        2 * 3.0/637.0,
                        3 * 3.0/637.0,
                        4 * 3.0/637.0,
                        5 * 3.0/637.0,
                        6 * 3.0/637.0,
                    });

                    DenseVector output = new DenseVector(6);
                    output.SetLength(3);
                    output.Assignment(input);

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }
            }

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

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

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

                    DenseVector output = input.Clone();

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(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, },
                        { 5, 10, 15, },
                        { 6, 12, 18, },
                    });

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

                    var expected = new DenseVector(new Scalar[] {
                        1 * 3.0/182,
                        2 * 3.0/182,
                        3 * 3.0/182,
                    });

                    DenseVector output = input.Clone();

                    var lls = new CompleteOrthogonalDecomposition(ref matrix);
                    lls.UnconstrainedLeastSquaresSolution(ref output);

                    CheckEqual(expected, output);
                }
            }
        }

        public class PreMatrixMultiply : UnitTestSharp.TestFixture
        {
            public void WrongSize()
            {
                var matrix = new RectangularMatrix(new Scalar[,] {
                    { -7 },
                });

                var decomp = new CompleteOrthogonalDecomposition(ref matrix);

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

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

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

                var decomp = new CompleteOrthogonalDecomposition(ref matrix);

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

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

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

                var decomp = new CompleteOrthogonalDecomposition(ref matrix);

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

                decomp.PremultiplyByMatrix(ref outputVector);

                CheckEqual(expected, outputVector);
            }

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

                var decomp = new CompleteOrthogonalDecomposition(ref matrix);

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

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

                var outputVector = new DenseVector(4);
                outputVector.SetLength(2);
                outputVector.Assignment(inputVector);

                decomp.PremultiplyByMatrix(ref outputVector);

                CheckEqual(expected, outputVector);
            }

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

                var decomp = new CompleteOrthogonalDecomposition(ref matrix);

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

                decomp.PremultiplyByMatrix(ref outputVector);

                CheckEqual(expected, outputVector);
            }
        }
    }
}