/// <summary>Performs the QR factorization.</summary>
    protected override void InternalCompute() 
    {
      int m = matrix.Rows;
      int n = matrix.Columns;
      
#if MANAGED
      int minmn = m < n ? m : n;
      r_ = new ComplexFloatMatrix(matrix); // create a copy
      ComplexFloatVector[] u = new ComplexFloatVector[minmn];
      for (int i = 0; i < minmn; i++) 
      {
        u[i] = Householder.GenerateColumn(r_, i, m - 1, i);
        Householder.UA(u[i], r_, i, m - 1, i + 1, n - 1);
      }
      q_ = ComplexFloatMatrix.CreateIdentity(m);
      for (int i = minmn - 1; i >= 0; i--) 
      {
        Householder.UA(u[i], q_, i, m - 1, i, m - 1);
      }
#else
      qr = ComplexFloatMatrix.ToLinearComplexArray(matrix);
      jpvt = new int[n];
      jpvt[0] = 1;
      Lapack.Geqp3.Compute(m, n, qr, m, jpvt, out tau);
      r_ = new ComplexFloatMatrix(m, n);
      // Populate R

      for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
          if (i <= j) {
            r_.data[j * m + i] = qr[(jpvt[j]-1) * m + i];
          }
          else {
            r_.data[j * m + i] = ComplexFloat.Zero;
          }
        }
      }

      q_ = new ComplexFloatMatrix(m, m);
      for (int i = 0; i < m; i++) {
        for (int j = 0; j < m; j++) {
          if (j < n)
            q_.data[j * m + i] = qr[j * m + i];
          else
            q_.data[j * m + i] = ComplexFloat.Zero;
        }
      }
      if( m < n ){
        Lapack.Ungqr.Compute(m, m, m, q_.data, m, tau);
      } else{
        Lapack.Ungqr.Compute(m, m, n, q_.data, m, tau);
      }
#endif
      for (int i = 0; i < m; i++) 
      {
        if (q_[i, i] == 0)
          isFullRank = false;
      }
    }

		public void CtorDimensions()
		{
			ComplexFloatVector test = new ComplexFloatVector(2);

			Assert.AreEqual(test.Length, 2);
			Assert.AreEqual(test[0], (ComplexFloat)0);
			Assert.AreEqual(test[1], (ComplexFloat)0);
		}

		public void CtorInitialValues()
		{
			ComplexFloatVector test = new ComplexFloatVector(2, (ComplexFloat)1);

			Assert.AreEqual(test.Length, 2);
			Assert.AreEqual(test[0], (ComplexFloat)1);
			Assert.AreEqual(test[1], (ComplexFloat)1);
		}

 public void CurrentException2() 
 {
   ComplexFloatVector test = new ComplexFloatVector(new ComplexFloat[2]{1f,2f});
   IEnumerator enumerator = test.GetEnumerator();
   enumerator.MoveNext();
   enumerator.MoveNext();
   enumerator.MoveNext();
   object value=enumerator.Current;
 }

		public void CtorArray()
		{
			float[] testvector = new float[2] { 0, 1 };

			ComplexFloatVector test = new ComplexFloatVector(testvector);
			Assert.AreEqual(test.Length, testvector.Length);
			Assert.AreEqual(test[0], (ComplexFloat)testvector[0]);
			Assert.AreEqual(test[1], (ComplexFloat)testvector[1]);
		}

 public void Current()
 {
   ComplexFloatVector test = new ComplexFloatVector(new ComplexFloat[2]{1f,2f});
   IEnumerator enumerator = test.GetEnumerator();
   bool movenextresult;
   
   movenextresult=enumerator.MoveNext();
   Assert.IsTrue(movenextresult);
   Assert.AreEqual(enumerator.Current,test[0]);
   
   movenextresult=enumerator.MoveNext();
   Assert.IsTrue(movenextresult);
   Assert.AreEqual(enumerator.Current,test[1]);
   
   movenextresult=enumerator.MoveNext();
   Assert.IsFalse(movenextresult);
 }

 public void ZeroLengthVectorTestsforConstructor1()
 {
   ComplexFloatVector cfv = new ComplexFloatVector(1);
   cfv.RemoveAt(0);
   ComplexFloatLevinson cfl = new ComplexFloatLevinson(cfv, cfv);
 }

 public void FirstElementTestforStaticInverse()
 {
   ComplexFloatVector cfv = new ComplexFloatVector(3, 1.0f);
   ComplexFloatMatrix X = ComplexFloatLevinson.Inverse(cfv, TR3);
 }

 public void SingularTestforStaticSolveMatrix()
 {
   ComplexFloatVector cfv = new ComplexFloatVector(3, 1.0f);
   ComplexFloatMatrix X = ComplexFloatLevinson.Solve(cfv, cfv, ComplexFloatMatrix.CreateIdentity(3));
 }

 public void ZeroVectorLengthTestforStaticSolveMatrix()
 {
   ComplexFloatVector LC = new ComplexFloatVector(1);
   LC.RemoveAt(0);
   ComplexFloatMatrix X = ComplexFloatLevinson.Solve(LC, TR10, ComplexFloatMatrix.CreateIdentity(10));
 }

 public void FirstElementTestforStaticSolveVector()
 {
   ComplexFloatVector cfv = new ComplexFloatVector(3, 1.0f);
   ComplexFloatVector X = ComplexFloatLevinson.Solve(cfv, TR3, Y3);
 }

    public void SetupTestCases()
    {
      // unit testing values - order 1

      LC1 = new ComplexFloatVector(1);
      LC1[0] = new ComplexFloat(+1.0000000E+000f, +1.0000000E+000f);

      TR1 = new ComplexFloatVector(1);
      TR1[0] = new ComplexFloat(+1.0000000E+000f, +1.0000000E+000f);

      L1 = new ComplexFloatMatrix(1);
      L1[0, 0] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);

      D1 = new ComplexFloatVector(1);
      D1[0] = new ComplexFloat(+5.0000000E-001f, -5.0000000E-001f);

      U1 = new ComplexFloatMatrix(1);
      U1[0, 0] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);

      Det1 = new ComplexFloat(+1.0000000E+000f, +1.0000000E+000f);

      I1 = new ComplexFloatMatrix(1);
      I1[0, 0] = new ComplexFloat(+5.0000000E-001f, -5.0000000E-001f);

      X1 = new ComplexFloatVector(1);
      X1[0] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);

      Y1 = new ComplexFloatVector(1);
      Y1[0] = new ComplexFloat(+1.0000000E+000f, +1.0000000E+000f);

      // unit testing values - order 2

      LC2 = new ComplexFloatVector(2);
      LC2[0] = new ComplexFloat(+3.0000000E+000f, +3.0000000E+000f);
      LC2[1] = new ComplexFloat(+2.0000000E+000f, +0.0000000E+000f);

      TR2 = new ComplexFloatVector(2);
      TR2[0] = new ComplexFloat(+3.0000000E+000f, +3.0000000E+000f);
      TR2[1] = new ComplexFloat(+2.0000000E+000f, +0.0000000E+000f);

      L2 = new ComplexFloatMatrix(2);
      L2[0, 0] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);
      L2[1, 0] = new ComplexFloat(-3.3333333E-001f, +3.3333333E-001f);
      L2[1, 1] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);

      D2 = new ComplexFloatVector(2);
      D2[0] = new ComplexFloat(+1.6666667E-001f, -1.6666667E-001f);
      D2[1] = new ComplexFloat(+1.2352941E-001f, -1.9411765E-001f);

      U2 = new ComplexFloatMatrix(2);
      U2[0, 0] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);
      U2[0, 1] = new ComplexFloat(-3.3333333E-001f, +3.3333333E-001f);
      U2[1, 1] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);

      Det2 = new ComplexFloat(-4.0000000E+000f, +1.8000000E+001f);

      I2 = new ComplexFloatMatrix(2);
      I2[0, 0] = new ComplexFloat(+1.2352941E-001f, -1.9411765E-001f);
      I2[0, 1] = new ComplexFloat(+2.3529412E-002f, +1.0588235E-001f);
      I2[1, 0] = new ComplexFloat(+2.3529412E-002f, +1.0588235E-001f);
      I2[1, 1] = new ComplexFloat(+1.2352941E-001f, -1.9411765E-001f);

      X2 = new ComplexFloatVector(2);
      X2[0] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);
      X2[1] = new ComplexFloat(+2.0000000E+000f, +0.0000000E+000f);

      Y2 = new ComplexFloatVector(2);
      Y2[0] = new ComplexFloat(+7.0000000E+000f, +3.0000000E+000f);
      Y2[1] = new ComplexFloat(+8.0000000E+000f, +6.0000000E+000f);

      // unit testing values - order 3

      LC3 = new ComplexFloatVector(3);
      LC3[0] = new ComplexFloat(+3.0000000E+000f, +3.0000000E+000f);
      LC3[1] = new ComplexFloat(+2.0000000E+000f, +0.0000000E+000f);
      LC3[2] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);

      TR3 = new ComplexFloatVector(3);
      TR3[0] = new ComplexFloat(+3.0000000E+000f, +3.0000000E+000f);
      TR3[1] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);
      TR3[2] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);

      L3 = new ComplexFloatMatrix(3);
      L3[0, 0] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);
      L3[1, 0] = new ComplexFloat(-3.3333333E-001f, +3.3333333E-001f);
      L3[1, 1] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);
      L3[2, 0] = new ComplexFloat(-1.7073171E-001f, -3.6585366E-002f);
      L3[2, 1] = new ComplexFloat(-2.9878049E-001f, +3.1097561E-001f);
      L3[2, 2] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);

      D3 = new ComplexFloatVector(3);
      D3[0] = new ComplexFloat(+1.6666667E-001f, -1.6666667E-001f);
      D3[1] = new ComplexFloat(+1.4634146E-001f, -1.8292683E-001f);
      D3[2] = new ComplexFloat(+1.4776571E-001f, -1.9120527E-001f);

      U3 = new ComplexFloatMatrix(3);
      U3[0, 0] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);
      U3[0, 1] = new ComplexFloat(-1.6666667E-001f, +1.6666667E-001f);
      U3[1, 1] = new ComplexFloat(+1.0000000E+000f, +0.0000000E+000f);
      U3[0, 2] = new ComplexFloat(-1.5243902E-001f, +1.2804878E-001f);
//.........这里部分代码省略.........

 public void StaticMultiplyMatrixNonConformVector()
 {
   ComplexFloatMatrix a = new ComplexFloatMatrix(2);
   ComplexFloatVector b = new ComplexFloatVector(3, 2.0f);
   ComplexFloatVector c = a * b;
 }

 public void StaticMultiplyMatrixVector()
 {
   ComplexFloatMatrix a = new ComplexFloatMatrix(2);
   a[0,0] = new ComplexFloat(1);
   a[0,1] = new ComplexFloat(2);
   a[1,0] = new ComplexFloat(3);
   a[1,1] = new ComplexFloat(4);
   ComplexFloatVector b = new ComplexFloatVector(2, 2.0f);
   ComplexFloatVector c = ComplexFloatMatrix.Multiply(a,b);
   Assert.AreEqual(c[0],new ComplexFloat(6));
   Assert.AreEqual(c[1],new ComplexFloat(14));
 }

 public void OperatorMultiplyNullMatrixVector()
 {
   ComplexFloatMatrix a = null;
   ComplexFloatVector b = new ComplexFloatVector(2, 2.0f);
   ComplexFloatVector c = a * b;
 }

    /// <overloads>
    /// Solve a symmetric square Toeplitz system.
    /// </overloads>
    /// <summary>
    /// Solve a symmetric square Toeplitz system with a right-side vector.
    /// </summary>
    /// <param name="T">The left-most column of the Toeplitz matrix.</param>
    /// <param name="Y">The right-side vector of the system.</param>
    /// <returns>The solution vector.</returns>
    /// <exception cref="ArgumentNullException">
    /// <B>T</B> and/or <B>Y</B> are null references
    /// </exception>
    /// <exception cref="RankException">
    /// The length of <B>T</B> does not match the length of <B>Y</B>.
    /// </exception>
    /// <exception cref="SingularMatrixException">
    /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
    /// </exception>
    /// <remarks>
    /// This method solves the linear system <B>AX</B> = <B>Y</B>. Where
    /// <B>T</B> is a symmetric square Toeplitz matrix, <B>X</B> is an unknown
    /// vector and <B>Y</B> is a known vector.
    /// <para>
    /// This static member combines the <b>UDL</b> decomposition and the calculation of the solution into a
    /// single algorithm. When compared to the non-static member it requires minimal data storage
    /// and suffers from no speed penalty.
    /// </para>
    /// </remarks>
    public static ComplexFloatVector Solve(IROComplexFloatVector T, IROComplexFloatVector Y)
    {

      ComplexFloatVector X;

      // check parameters
      if (T == null)
      {
        throw new System.ArgumentNullException("T");
      }
      else if (Y == null)
      {
        throw new System.ArgumentNullException("Y");
      }
      else if (T.Length != Y.Length)
      {
        throw new RankException("The length of T and Y are not equal.");
      }
      else
      {

        // allocate memory
        int N = T.Length;
        X = new ComplexFloatVector(N);                    // solution vector
        ComplexFloat e;                                   // prediction error

        // setup zero order solution
        e = T[0];
        if (e == ComplexFloat.Zero)
        {
          throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
        }
        X[0] = Y[0] / T[0];

        if (N > 1)
        {
          ComplexFloatVector a = new ComplexFloatVector(N - 1);   // prediction coefficients
          ComplexFloatVector Z = new ComplexFloatVector(N - 1);   // temporary storage vector
          ComplexFloat g;                                   // reflection coefficient
          ComplexFloat inner;                               // inner product
          ComplexFloat k;
          int i, j, l;

          // calculate solution for successive orders
          for (i = 1; i < N; i++)
          {

            // calculate first inner product
            inner = T[i];
            for (j = 0, l = i - 1; j < i - 1; j++, l--)
            {
              inner += a[j] * T[l];
            }

            // update predictor coefficients
            g = -(inner / e);
            for (j = 0, l = i - 2; j < i - 1; j++, l--)
            {
              Z[j] = a[j] + g * a[l];
            }

            // copy vector
            for (j = 0; j < i - 1; j++)
            {
              a[j] = Z[j];
            }

            a[i - 1] = g;

            e *= (ComplexFloat.One - g * g);
            if (e == ComplexFloat.Zero)
            {
//.........这里部分代码省略.........

    public void SingularityPropertyTest2()
    {
      ComplexFloatVector LC = new ComplexFloatVector(4);
      LC[0] = new ComplexFloat(4.0f);
      LC[1] = new ComplexFloat(2.0f);
      LC[2] = new ComplexFloat(1.0f);
      LC[3] = new ComplexFloat(0.0f);

      ComplexFloatVector TR = new ComplexFloatVector(4);
      TR[0] = new ComplexFloat(4.0f);
      TR[1] = new ComplexFloat(8.0f);
      TR[2] = new ComplexFloat(2.0f);
      TR[3] = new ComplexFloat(1.0f);

      ComplexFloatLevinson cfl = new ComplexFloatLevinson(LC,TR);
      Assert.IsTrue(cfl.IsSingular);
    }

 public void StaticMultiplyNullMatrixVector()
 {
   ComplexFloatMatrix a = null;
   ComplexFloatVector b = new ComplexFloatVector(2, 2.0f);
   ComplexFloatVector c = ComplexFloatMatrix.Multiply(a,b);
 }

 public void ZeroVectorLengthTestforStaticSolveVector()
 {
   ComplexFloatVector LC = new ComplexFloatVector(1);
   LC.RemoveAt(0);
   ComplexFloatVector X = ComplexFloatLevinson.Solve(LC, TR10, Y10);
 }

 public void MemberMultiplyVector()
 {
   ComplexFloatMatrix a = new ComplexFloatMatrix(2);
   a[0,0] = new ComplexFloat(1);
   a[0,1] = new ComplexFloat(2);
   a[1,0] = new ComplexFloat(3);
   a[1,1] = new ComplexFloat(4);
   ComplexFloatVector b = new ComplexFloatVector(2, 2.0f);
   a.Multiply(b);
   Assert.AreEqual(a[0,0],new ComplexFloat(6));
   Assert.AreEqual(a[1,0],new ComplexFloat(14));
   Assert.AreEqual(a.ColumnLength, 1);
   Assert.AreEqual(a.RowLength, 2);
 }