public void CanCreateSparseVectorFromArray()
        {
            var data = new Complex32[Data.Length];
            Array.Copy(Data, data, Data.Length);
            var vector = new SparseVector(data);

            for (var i = 0; i < data.Length; i++)
            {
                Assert.AreEqual(data[i], vector[i]);
            }
        }

        public void CanWriteMatrices()
        {
            Matrix mat1 = new DenseMatrix(5, 3);
            for (var i = 0; i < mat1.ColumnCount; i++)
            {
                mat1[i, i] = new Complex32(i + .1f, i + .1f);
            }

            Matrix mat2 = new DenseMatrix(4, 5);
            for (var i = 0; i < mat2.RowCount; i++)
            {
                mat2[i, i] = new Complex32(i + .1f, i + .1f);
            }

            Matrix mat3 = new SparseMatrix(5, 4);
            for (var i = 0; i < mat3.ColumnCount; i++)
            {
                mat3[i, i] = new Complex32(i + .1f, i + .1f);
            }

            Matrix mat4 = new SparseMatrix(3, 5);
            for (var i = 0; i < mat4.RowCount; i++)
            {
                mat4[i, i] = new Complex32(i + .1f, i + .1f);
            }

            var write = new[] { mat1, mat2, mat3, mat4 };

            var names = new[] { "mat1", "dense_matrix_2", "s1", "sparse2" };
            if (File.Exists("test.mat"))
            {
                File.Delete("test.mat");
            }

            var writer = new MatlabMatrixWriter("test.mat");
            writer.WriteMatrices(write, names);
            writer.Dispose();

            var reader = new MatlabMatrixReader("test.mat");
            var read = reader.ReadMatrices(names);

            Assert.AreEqual(write.Length, read.Count);

            for (var i = 0; i < write.Length; i++)
            {
                var w = write[i];
                var r = read[names[i]];

                Assert.AreEqual(w.RowCount, r.RowCount);
                Assert.AreEqual(w.ColumnCount, r.ColumnCount);
                Assert.IsTrue(w.Equals(r));
            }
        }

        public void CanCreateDenseVectorFromArray()
        {
            var data = new Complex32[Data.Length];
            Array.Copy(Data, data, Data.Length);
            var vector = new DenseVector(data);

            for (var i = 0; i < data.Length; i++)
            {
                Assert.AreEqual(data[i], vector[i]);
            }

            vector[0] = new Complex32(10.0f, 1);
            Assert.AreEqual(new Complex32(10.0f, 1), data[0]);
        }

        public void CanMultiplyWithComplex([Values(0, 1, 2.2f)] float real)
        {
            var value = new Complex32(real, 1.0f);
            var matrix = TestMatrices["Singular3x3"];
            var clone = matrix.Clone();
            clone = clone.Multiply(value);

            for (var i = 0; i < matrix.RowCount; i++)
            {
                for (var j = 0; j < matrix.ColumnCount; j++)
                {
                    Assert.AreEqual(matrix[i, j] * value, clone[i, j]);
                }
            }
        }

        public void CanSetDiagonalVector([Values("Square3x3", "Wide2x3", "Tall3x2")] string name, [Values(new float[] { 1, 2, 3 }, new float[] { 1, 2 }, new float[] { 1, 2 })] float[] real)
        {
            var diagonal = new Complex32[real.Length];
            for (var i = 0; i < real.Length; i++)
            {
                diagonal[i] = new Complex32(real[i], 1);
            }

            var matrix = TestMatrices[name];
            var vector = CreateVector(diagonal);
            matrix.SetDiagonal(vector);

            var min = Math.Min(matrix.ColumnCount, matrix.RowCount);
            Assert.AreEqual(diagonal.Length, min);

            for (var i = 0; i < vector.Count; i++)
            {
                Assert.AreEqual(vector[i], matrix[i, i]);
            }
        }

        public void CheckSparseMechanismBySettingValues()
        {
            var vector = new SparseVector(10000);

            // Add non-zero elements
            vector[200] = new Complex32(1.5f, 1);
            Assert.AreEqual(new Complex32(1.5f, 1), vector[200]);
            Assert.AreEqual(1, vector.NonZerosCount);

            vector[500] = new Complex32(3.5f, 1);
            Assert.AreEqual(new Complex32(3.5f, 1), vector[500]);
            Assert.AreEqual(2, vector.NonZerosCount);

            vector[800] = new Complex32(5.5f, 1);
            Assert.AreEqual(new Complex32(5.5f, 1), vector[800]);
            Assert.AreEqual(3, vector.NonZerosCount);

            vector[0] = new Complex32(7.5f, 1);
            Assert.AreEqual(new Complex32(7.5f, 1), vector[0]);
            Assert.AreEqual(4, vector.NonZerosCount);

            // Remove non-zero elements
            vector[200] = Complex32.Zero;
            Assert.AreEqual(Complex32.Zero, vector[200]);
            Assert.AreEqual(3, vector.NonZerosCount);

            vector[500] = Complex32.Zero;
            Assert.AreEqual(Complex32.Zero, vector[500]);
            Assert.AreEqual(2, vector.NonZerosCount);

            vector[800] = Complex32.Zero;
            Assert.AreEqual(Complex32.Zero, vector[800]);
            Assert.AreEqual(1, vector.NonZerosCount);

            vector[0] = Complex32.Zero;
            Assert.AreEqual(Complex32.Zero, vector[0]);
            Assert.AreEqual(0, vector.NonZerosCount);
        }

        public void CheckSparseMechanismByZeroMultiply()
        {
            var vector = new SparseVector(10000);

            // Add non-zero elements
            vector[200] = new Complex32(1.5f, 1);
            vector[500] = new Complex32(3.5f, 1);
            vector[800] = new Complex32(5.5f, 1);
            vector[0] = new Complex32(7.5f, 1);

            // Multiply by 0
            vector *= 0;
            Assert.AreEqual(Complex32.Zero, vector[200]);
            Assert.AreEqual(Complex32.Zero, vector[500]);
            Assert.AreEqual(Complex32.Zero, vector[800]);
            Assert.AreEqual(Complex32.Zero, vector[0]);
            Assert.AreEqual(0, vector.NonZerosCount);
        }

        public void CanSetDiagonalVector(string name, float[] real)
        {
            var diagonal = new Complex32[real.Length];
            for (var i = 0; i < real.Length; i++)
            {
                diagonal[i] = new Complex32(real[i], 1);
            }

            var matrix = TestMatrices[name];
            var vector = CreateVector(diagonal);
            matrix.SetDiagonal(vector);

            var min = Math.Min(matrix.ColumnCount, matrix.RowCount);
            Assert.AreEqual(diagonal.Length, min);

            for (var i = 0; i < vector.Count; i++)
            {
                Assert.AreEqual(vector[i], matrix[i, i]);
            }
        }

        public void CanMultiplySparseVectorByComplexUsingOperators()
        {
            var vector = new SparseVector(Data);
            vector = vector * new Complex32(2.0f, 1);

            for (var i = 0; i < Data.Length; i++)
            {
                Assert.AreEqual(Data[i] * new Complex32(2.0f, 1), vector[i]);
            }

            vector = vector * 1.0f;
            for (var i = 0; i < Data.Length; i++)
            {
                Assert.AreEqual(Data[i] * new Complex32(2.0f, 1), vector[i]);
            }

            vector = new SparseVector(Data);
            vector = new Complex32(2.0f, 1) * vector;

            for (var i = 0; i < Data.Length; i++)
            {
                Assert.AreEqual(Data[i] * new Complex32(2.0f, 1), vector[i]);
            }

            vector = 1.0f * vector;
            for (var i = 0; i < Data.Length; i++)
            {
                Assert.AreEqual(Data[i] * new Complex32(2.0f, 1), vector[i]);
            }
        }

 /// <summary>Sets the <paramref name="value"/> at the given <paramref name="index"/>.</summary>
 /// <param name="index">The index of the value to get or set.</param>
 /// <param name="value">The value to set.</param>
 protected internal override void At(int index, Complex32 value)
 {
     this[index] = value;
 }

 public void CanSum()
 {
     Complex32[] testData = { new Complex32(-20, -1), new Complex32(-10, -1), new Complex32(10, 1), new Complex32(20, 1), new Complex32(30, 1) };
     var vector = CreateVector(testData);
     var actual = vector.Sum();
     var expected = new Complex32(30, 1);
     Assert.AreEqual(expected, actual);
 }

        /// <summary>
        /// Nonsymmetric reduction from Hessenberg to real Schur form.
        /// </summary>
        /// <param name="vectorV">Data array of the eigenvectors</param>
        /// <param name="dataEv">Data array of matrix V (eigenvectors)</param>
        /// <param name="matrixH">Array for internal storage of nonsymmetric Hessenberg form.</param>
        /// <param name="order">Order of initial matrix</param>
        /// <remarks>This is derived from the Algol procedure hqr2,
        /// by Martin and Wilkinson, Handbook for Auto. Comp.,
        /// Vol.ii-Linear Algebra, and the corresponding
        /// Fortran subroutine in EISPACK.</remarks>
        internal static void NonsymmetricReduceHessenberToRealSchur(Numerics.Complex32[] vectorV, Numerics.Complex32[] dataEv, Numerics.Complex32[] matrixH, int order)
        {
            // Initialize
            var n = order - 1;
            var eps = (float) Precision.SingleMachinePrecision;

            float norm;
            Numerics.Complex32 x, y, z, exshift = Numerics.Complex32.Zero;

            // Outer loop over eigenvalue index
            var iter = 0;
            while (n >= 0)
            {
                // Look for single small sub-diagonal element
                var l = n;
                while (l > 0)
                {
                    var lm1 = l - 1;
                    var lm1O = lm1*order;
                    var lO = l*order;
                    var tst1 = Math.Abs(matrixH[lm1O + lm1].Real) + Math.Abs(matrixH[lm1O + lm1].Imaginary) + Math.Abs(matrixH[lO + l].Real) + Math.Abs(matrixH[lO + l].Imaginary);
                    if (Math.Abs(matrixH[lm1O + l].Real) < eps*tst1)
                    {
                        break;
                    }

                    l--;
                }

                var nm1 = n - 1;
                var nm1O = nm1*order;
                var nO = n*order;
                var nOn = nO + n;
                // Check for convergence
                // One root found
                if (l == n)
                {
                    matrixH[nOn] += exshift;
                    vectorV[n] = matrixH[nOn];
                    n--;
                    iter = 0;
                }
                else
                {
                    // Form shift
                    Numerics.Complex32 s;
                    if (iter != 10 && iter != 20)
                    {
                        s = matrixH[nOn];
                        x = matrixH[nO + nm1]*matrixH[nm1O + n].Real;

                        if (x.Real != 0.0f || x.Imaginary != 0.0f)
                        {
                            y = (matrixH[nm1O + nm1] - s)/2.0f;
                            z = ((y*y) + x).SquareRoot();
                            if ((y.Real*z.Real) + (y.Imaginary*z.Imaginary) < 0.0)
                            {
                                z *= -1.0f;
                            }

                            x /= y + z;
                            s = s - x;
                        }
                    }
                    else
                    {
                        // Form exceptional shift
                        s = Math.Abs(matrixH[nm1O + n].Real) + Math.Abs(matrixH[(n - 2)*order + nm1].Real);
                    }

                    for (var i = 0; i <= n; i++)
                    {
                        matrixH[i*order + i] -= s;
                    }

                    exshift += s;
                    iter++;

                    // Reduce to triangle (rows)
                    for (var i = l + 1; i <= n; i++)
                    {
                        var im1 = i - 1;
                        var im1O = im1*order;
                        var im1Oim1 = im1O + im1;
                        s = matrixH[im1O + i].Real;
                        norm = SpecialFunctions.Hypotenuse(matrixH[im1Oim1].Magnitude, s.Real);
                        x = matrixH[im1Oim1]/norm;
                        vectorV[i - 1] = x;
                        matrixH[im1Oim1] = norm;
//.........这里部分代码省略.........

        /// <summary>
        /// Solves a system of linear equations, <b>Ax = b</b>, with A EVD factorized.
        /// </summary>
        /// <param name="input">The right hand side vector, <b>b</b>.</param>
        /// <param name="result">The left hand side <see cref="Matrix{T}"/>, <b>x</b>.</param>
        public override void Solve(Vector<Numerics.Complex32> input, Vector<Numerics.Complex32> result)
        {
            if (input == null)
            {
                throw new ArgumentNullException("input");
            }

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

            // Ax=b where A is an m x m matrix
            // Check that b is a column vector with m entries
            if (VectorEv.Count != input.Count)
            {
                throw new ArgumentException(Resources.ArgumentVectorsSameLength);
            }

            // Check that x is a column vector with n entries
            if (VectorEv.Count != result.Count)
            {
                throw new ArgumentException(Resources.ArgumentMatrixDimensions);
            }

            if (IsSymmetric)
            {
                // Symmetric case -> x = V * inv(λ) * VH * b;
                var order = VectorEv.Count;
                var tmp = new Numerics.Complex32[order];
                Numerics.Complex32 value;

                for (var j = 0; j < order; j++)
                {
                    value = 0;
                    if (j < order)
                    {
                        for (var i = 0; i < order; i++)
                        {
                            value += ((DenseMatrix) MatrixEv).Values[(j*order) + i].Conjugate()*input[i];
                        }

                        value /= (float) VectorEv[j].Real;
                    }

                    tmp[j] = value;
                }

                for (var j = 0; j < order; j++)
                {
                    value = 0;
                    for (var i = 0; i < order; i++)
                    {
                        value += ((DenseMatrix) MatrixEv).Values[(i*order) + j]*tmp[i];
                    }

                    result[j] = value;
                }
            }
            else
            {
                throw new ArgumentException(Resources.ArgumentMatrixSymmetric);
            }
        }

        public virtual void CanSetColumnWithArray([Values("Singular3x3", "Square3x3", "Tall3x2", "Wide2x3")] string name, [Values(new float[] { 1, 2, 3 }, new float[] { 1, 2, 3 }, new float[] { 1, 2, 3 }, new float[] { 1, 2 })] float[] real)
        {
            var column = new Complex32[real.Length];
            for (var i = 0; i < real.Length; i++)
            {
                column[i] = new Complex32(real[i], 1);
            }

            var matrix = TestMatrices[name];
            for (var i = 0; i < matrix.ColumnCount; i++)
            {
                matrix.SetColumn(i, column);
                for (var j = 0; j < matrix.RowCount; j++)
                {
                    Assert.AreEqual(matrix[j, i], column[j]);
                }
            }
        }

        /// <summary>
        /// Reduces a complex hermitian matrix to a real symmetric tridiagonal matrix using unitary similarity transformations.
        /// </summary>
        /// <param name="matrixA">Source matrix to reduce</param>
        /// <param name="d">Output: Arrays for internal storage of real parts of eigenvalues</param>
        /// <param name="e">Output: Arrays for internal storage of imaginary parts of eigenvalues</param>
        /// <param name="tau">Output: Arrays that contains further information about the transformations.</param>
        /// <param name="order">Order of initial matrix</param>
        /// <remarks>This is derived from the Algol procedures HTRIDI by 
        /// Smith, Boyle, Dongarra, Garbow, Ikebe, Klema, Moler, and Wilkinson, Handbook for 
        /// Auto. Comp., Vol.ii-Linear Algebra, and the corresponding 
        /// Fortran subroutine in EISPACK.</remarks>
        internal static void SymmetricTridiagonalize(Numerics.Complex32[] matrixA, float[] d, float[] e, Numerics.Complex32[] tau, int order)
        {
            float hh;
            tau[order - 1] = Numerics.Complex32.One;

            for (var i = 0; i < order; i++)
            {
                d[i] = matrixA[i*order + i].Real;
            }

            // Householder reduction to tridiagonal form.
            for (var i = order - 1; i > 0; i--)
            {
                // Scale to avoid under/overflow.
                var scale = 0.0f;
                var h = 0.0f;

                for (var k = 0; k < i; k++)
                {
                    scale = scale + Math.Abs(matrixA[k*order + i].Real) + Math.Abs(matrixA[k*order + i].Imaginary);
                }

                if (scale == 0.0f)
                {
                    tau[i - 1] = Numerics.Complex32.One;
                    e[i] = 0.0f;
                }
                else
                {
                    for (var k = 0; k < i; k++)
                    {
                        matrixA[k*order + i] /= scale;
                        h += matrixA[k*order + i].MagnitudeSquared;
                    }

                    Numerics.Complex32 g = (float) Math.Sqrt(h);
                    e[i] = scale*g.Real;

                    Numerics.Complex32 temp;
                    var im1Oi = (i - 1)*order + i;
                    var f = matrixA[im1Oi];
                    if (f.Magnitude != 0.0f)
                    {
                        temp = -(matrixA[im1Oi].Conjugate()*tau[i].Conjugate())/f.Magnitude;
                        h += f.Magnitude*g.Real;
                        g = 1.0f + (g/f.Magnitude);
                        matrixA[im1Oi] *= g;
                    }
                    else
                    {
                        temp = -tau[i].Conjugate();
                        matrixA[im1Oi] = g;
                    }

                    if ((f.Magnitude == 0.0f) || (i != 1))
                    {
                        f = Numerics.Complex32.Zero;
                        for (var j = 0; j < i; j++)
                        {
                            var tmp = Numerics.Complex32.Zero;
                            var jO = j*order;
                            // Form element of A*U.
                            for (var k = 0; k <= j; k++)
                            {
                                tmp += matrixA[k*order + j]*matrixA[k*order + i].Conjugate();
                            }

                            for (var k = j + 1; k <= i - 1; k++)
                            {
                                tmp += matrixA[jO + k].Conjugate()*matrixA[k*order + i].Conjugate();
                            }

                            // Form element of P
                            tau[j] = tmp/h;
                            f += (tmp/h)*matrixA[jO + i];
                        }

                        hh = f.Real/(h + h);

                        // Form the reduced A.
                        for (var j = 0; j < i; j++)
                        {
                            f = matrixA[j*order + i].Conjugate();
                            g = tau[j] - (hh*f);
                            tau[j] = g.Conjugate();

                            for (var k = 0; k <= j; k++)
                            {
//.........这里部分代码省略.........

 public void SetDiagonalArrayWithUnequalLengthThrowsArgumentException()
 {
     var array = new Complex32[] { 1, 2, 3 };
     Assert.Throws<ArgumentException>(() => TestMatrices["Wide2x3"].SetDiagonal(array));
 }

 public virtual void SetColumnWithArrayUnequalLengthThrowsArgumentException()
 {
     var array = new Complex32[] { 1, 2, 3, 4, 5 };
     Assert.Throws<ArgumentException>(() => TestMatrices["Singular3x3"].SetColumn(1, array));
 }

 /// <summary>
 /// Creates a vector from an array.
 /// </summary>
 /// <param name="data">The array to create this vector from.</param>
 /// <returns>The new vector. </returns>
 protected override Vector CreateVector(Complex32[] data)
 {
     return SparseVector.OfEnumerable(data);
 }

 /// <summary>
 /// Creates a matrix from a 2D array.
 /// </summary>
 /// <param name="data">The 2D array to create this matrix from.</param>
 /// <returns>A matrix with the given values.</returns>
 protected override Matrix CreateMatrix(Complex32[,] data)
 {
     return SparseMatrix.OfArray(data);
 }

        public virtual void CanSetRowWithVector(string name, float[] real)
        {
            var row = new Complex32[real.Length];
            for (var i = 0; i < real.Length; i++)
            {
                row[i] = new Complex32(real[i], 1);
            }

            var matrix = TestMatrices[name];
            var rowVector = CreateVector(row);
            for (var i = 0; i < matrix.RowCount; i++)
            {
                matrix.SetRow(i, row);
                for (var j = 0; j < matrix.ColumnCount; j++)
                {
                    Assert.AreEqual(matrix[i, j], rowVector[j]);
                }
            }
        }