public void CtorDimensions()
 {
   ComplexDoubleVector test = new ComplexDoubleVector(2);
   
   Assert.AreEqual(test.Length, 2);
   Assert.AreEqual(test[0],new Complex(0));
   Assert.AreEqual(test[1],new Complex(0));
 }

 public void CtorInitialValues()
 {
   ComplexDoubleVector test = new ComplexDoubleVector(2, new Complex(1,-1));
   
   Assert.AreEqual(test.Length, 2);
   Assert.AreEqual(test[0],new Complex(1,-1));
   Assert.AreEqual(test[1],new Complex(1,-1));
 }

    /// <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 ComplexDoubleMatrix(matrix); // create a copy
      ComplexDoubleVector[] u = new ComplexDoubleVector[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_ = ComplexDoubleMatrix.CreateIdentity(m);
      for (int i = minmn - 1; i >= 0; i--) 
      {
        Householder.UA(u[i], q_, i, m - 1, i, m - 1);
      }
#else
      qr = ComplexDoubleMatrix.ToLinearComplexArray(matrix);
      jpvt = new int[n];
      jpvt[0] = 1;
      Lapack.Geqp3.Compute(m, n, qr, m, jpvt, out tau);
      r_ = new ComplexDoubleMatrix(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] = Complex.Zero;
          }
        }
      }

      q_ = new ComplexDoubleMatrix(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] = Complex.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 CurrentException2()
		{
			ComplexDoubleVector test = new ComplexDoubleVector(new Complex[2] { 1, 2 });
			IEnumerator enumerator = test.GetEnumerator();
			enumerator.MoveNext();
			enumerator.MoveNext();
			enumerator.MoveNext();
			object value = enumerator.Current;
		}

 public void CtorArray()
 {
   double[] testvector = new double[2]{0,1};
   
   ComplexDoubleVector test = new ComplexDoubleVector(testvector);
   Assert.AreEqual(test.Length,testvector.Length);
   Assert.AreEqual(test[0],new Complex(testvector[0]));
   Assert.AreEqual(test[1],new Complex(testvector[1]));
 }

		public static ComplexDoubleVector GenerateColumn(IComplexDoubleMatrix A, int r1, int r2, int c)
		{
			int ru = r2 - r1 + 1;
			ComplexDoubleVector u = new ComplexDoubleVector(r2 - r1 + 1);

			for (int i = r1; i <= r2; i++)
			{
				u[i - r1] = A[i, c];
				A[i, c] = Complex.Zero;
			}

			double norm = u.GetNorm();

			if (r1 == r2 || norm == 0)
			{
				A[r1, c] = new Complex(-u[0]);
				u[0] = System.Math.Sqrt(2);
				return u;
			}

			Complex scale = new Complex(1 / norm, 0);

			Complex t = Complex.Zero;
			Complex t1 = Complex.Zero;
			if (u[0].Real != 0 || u[0].Imag != 0)
			{
				t = u[0];
				t1 = ComplexMath.Conjugate(u[0]);
				t = ComplexMath.Absolute(t);
				t = t1 / t;
				scale = scale * t;
			}
			A[r1, c] = -Complex.One / scale;

			for (int i = 0; i < ru; i++)
			{
				u[i] = u[i] * scale;
			}

			u[0] = new Complex(u[0].Real + 1, 0);
			double s = System.Math.Sqrt(1 / u[0].Real);

			for (int i = 0; i < ru; i++)
			{
				u[i] = new Complex(s * u[i].Real, s * u[i].Imag);
			}
			return u;
		}

		public void Current()
		{
			ComplexDoubleVector test = new ComplexDoubleVector(new Complex[2] { 1, 2 });
			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 ICollection()
 {
   ComplexDoubleVector a = new ComplexDoubleVector(new Complex[4]{0,1,2,3});
   Complex[] b = new Complex[5];
   
   Assert.AreEqual(a.Count,a.Length);
   
   a.CopyTo(b,1);
   Assert.AreEqual(b[0],(Complex)0);
   Assert.AreEqual(b[1],(Complex)0); 
   Assert.AreEqual(b[2],(Complex)1);
   Assert.AreEqual(b[3],(Complex)2);
   Assert.AreEqual(b[4],(Complex)3); 
 }

 public void GetEnumerator()
 {
   ComplexDoubleVector a = new ComplexDoubleVector(new Complex[4]{0,1,2,3});
   IEnumerator dve = a.GetEnumerator();
   Complex b;
   bool c;
   
   c = dve.MoveNext();
   b = (Complex)dve.Current;
   Assert.AreEqual(c,true);
   Assert.AreEqual(b,(Complex)0);
   
   c = dve.MoveNext();
   b = (Complex)dve.Current;
   Assert.AreEqual(c,true);
   Assert.AreEqual(b,(Complex)1);
   
   c = dve.MoveNext();
   b = (Complex)dve.Current;
   Assert.AreEqual(c,true);
   Assert.AreEqual(b,(Complex)2);
   
   c = dve.MoveNext();
   b = (Complex)dve.Current;
   Assert.AreEqual(c,true);
   Assert.AreEqual(b,(Complex)3);
   
   c = dve.MoveNext();
   Assert.AreEqual(c,false);
 }

 public void Clone()
 {
   ComplexDoubleVector a = new ComplexDoubleVector(new double[4]{0,1,2,3});
   ComplexDoubleVector b = a.Clone();
   
   Assert.AreEqual(a[0],b[0]);
   Assert.AreEqual(a[1],b[1]);
   Assert.AreEqual(a[2],b[2]);
   Assert.AreEqual(a[3],b[3]);
   
   a=a*2;
   
   Assert.AreEqual(a[0],b[0]*2);
   Assert.AreEqual(a[1],b[1]*2);
   Assert.AreEqual(a[2],b[2]*2);
   Assert.AreEqual(a[3],b[3]*2);
 }

 public void GetEnumeratorException()
 {
   ComplexDoubleVector a = new ComplexDoubleVector(new Complex[4]{0,1,2,3});
   IEnumerator dve = a.GetEnumerator();
   
   Complex b = (Complex)dve.Current;
 }

 public void GetSum()
 {
   ComplexDoubleVector a = new ComplexDoubleVector(new double[4]{0,1,2,3});
   ComplexDoubleVector b = new ComplexDoubleVector(new double[4]{4,5,6,7});
   
   Assert.AreEqual(a.GetSum(),(Complex)6);
   Assert.AreEqual(a.GetSumMagnitudes(),6);
   
   Assert.AreEqual(b.GetSum(),(Complex)22);
   Assert.AreEqual(b.GetSumMagnitudes(),22);
 }

 public void Divide()
 {
   ComplexDoubleVector a = new ComplexDoubleVector(new double[4]{0,1,2,3});
   ComplexDoubleVector c = new ComplexDoubleVector(a);
   ComplexDoubleVector d = new ComplexDoubleVector(a);
   double scal = -4;
     
   c = a/scal;
   d = ComplexDoubleVector.Divide(a,scal);
     
   Assert.AreEqual(c[0],a[0]/scal);
   Assert.AreEqual(c[1],a[1]/scal);
   Assert.AreEqual(c[2],a[2]/scal);
   Assert.AreEqual(c[3],a[3]/scal);  
   
   Assert.AreEqual(d[0],a[0]/scal);
   Assert.AreEqual(d[1],a[1]/scal);
   Assert.AreEqual(d[2],a[2]/scal);
   Assert.AreEqual(d[3],a[3]/scal);
 }

    /// <summary>
    /// Solve the Yule-Walker equations for a symmetric square Toeplitz system
    /// </summary>
    /// <param name="R">The left-most column of the Toeplitz matrix.</param>
    /// <returns>The solution vector.</returns>
    /// <exception cref="ArgumentNullException">
    /// <B>R</B> is a null reference.
    /// </exception>
    /// <exception cref="ArgumentOutOfRangeException">
    /// The length of <B>R</B> must be greater than one.
    /// </exception>
    /// <exception cref="SingularMatrixException">
    /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
    /// </exception>
    /// <remarks>
    /// This member is used to solve the Yule-Walker system <B>AX</B> = -<B>a</B>,
    /// where <B>A</B> is a symmetric square Toeplitz matrix, constructed
    /// from the elements <B>R[0]</B>, ..., <B>R[N-2]</B> and
    /// the vector <B>a</B> is constructed from the elements
    /// <B>R[1]</B>, ..., <B>R[N-1]</B>.
    /// <para>
    /// Durbin's algorithm is used to solve the linear system. It requires
    /// approximately the <b>N</b> squared FLOPS to calculate the
    /// solution (<b>N</b> is the matrix order).
    /// </para>
    /// </remarks>
    public static ComplexDoubleVector YuleWalker(IROComplexDoubleVector R)
    {

      ComplexDoubleVector a;

      // check parameters
      if (R == null)
      {
        throw new System.ArgumentNullException("R");
      }
      else if (R.Length < 2)
      {
        throw new System.ArgumentOutOfRangeException("R", "The length of R must be greater than 1.");
      }
      else
      {

        int N = R.Length - 1;
        a = new ComplexDoubleVector(N);             // prediction coefficients
        ComplexDoubleVector Z = new ComplexDoubleVector(N);   // temporary storage vector
        Complex e;                    // predictor error
        Complex inner;                  // inner product
        Complex g;                    // reflection coefficient
        int i, j, l;

        // setup first order solution
        e = R[0];
        if (e == Complex.Zero)
        {
          throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
        }
        g = -R[1] / R[0];
        a[0] = g;

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

          e *= (Complex.One - g * g);
          if (e == Complex.Zero)
          {
            throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
          }

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

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

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

          a[i] = g;

        }
      }

      return a;
    }

    /// <overloads>
    /// Solve a symmetric square Toeplitz system.
    /// </overloads>
    /// <summary>
    /// Solve a symmetric square Toeplitz system with a right-side vector.
    /// </summary>
    /// <param name="Y">The right-hand side of the system.</param>
    /// <returns>The solution vector.</returns>
    /// <exception cref="ArgumentNullException">
    /// Parameter <B>Y</B> is a null reference.
    /// </exception>
    /// <exception cref="RankException">
    /// The length of <B>Y</B> is not equal to the number of rows in the Toeplitz matrix.
    /// </exception>
    /// <exception cref="SingularMatrixException">
    /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
    /// </exception>
    /// <remarks>
    /// This member solves the linear system <B>TX</B> = <B>Y</B>, where <B>T</B> is
    /// the symmetric square Toeplitz matrix, <B>X</B> is the unknown solution vector
    /// and <B>Y</B> is a known vector.
    /// <para>
    /// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
    /// using the Levinson algorithm, and then calculates the solution vector.
    /// </para>
    /// </remarks>
    public ComplexDoubleVector Solve(IROComplexDoubleVector Y)
    {
      ComplexDoubleVector X;

      // check parameters
      if (Y == null)
      {
        throw new System.ArgumentNullException("Y");
      }
      else if (m_Order != Y.Length)
      {
        throw new RankException("The length of Y is not equal to the number of rows in the Toeplitz matrix.");
      }

      Compute();

      if (m_IsSingular == true)
      {
        throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
      }

      int i, j, l;      // index/loop variables
      Complex Inner;      // inner product
      Complex G;        // scaling constant
      Complex[] A;        // reference to current order coefficients

      // allocate memory for solution
      X = new ComplexDoubleVector(m_Order);

      // setup zero order solution
      X[0] = Y[0] / m_LeftColumn[0];

      // solve systems of increasing order
      for (i = 1; i < m_Order; i++)
      {
        // calculate inner product
        Inner = Y[i];
        for (j = 0, l = i; j < i; j++, l--)
        {
          Inner -= X[j] * m_LeftColumn[l];
        }
        // get the current predictor coefficients row
        A = m_LowerTriangle[i];

        // update the solution vector
        G = Inner * m_Diagonal[i];
        for (j = 0; j <= i; j++)
        {
          X[j] += G * A[j];
        }

      }

      return X;
    }

 public void Add()
 {
   ComplexDoubleVector a = new ComplexDoubleVector(new double[4]{0,1,2,3});
   ComplexDoubleVector b = new ComplexDoubleVector(new double[4]{4,5,6,7});
   ComplexDoubleVector c = new ComplexDoubleVector(a.Length);
   ComplexDoubleVector d = new ComplexDoubleVector(b.Length);
   
   c = a+b;
   d = ComplexDoubleVector.Add(a,b);
   
   Assert.AreEqual(c[0],a[0]+b[0]);
   Assert.AreEqual(c[1],a[1]+b[1]);
   Assert.AreEqual(c[2],a[2]+b[2]);
   Assert.AreEqual(c[3],a[3]+b[3]);
   
   Assert.AreEqual(d[0],c[0]);
   Assert.AreEqual(d[1],c[1]);
   Assert.AreEqual(d[2],c[2]);
   Assert.AreEqual(d[3],c[3]);
   
   a.Add(b);
   
   Assert.AreEqual(c[0],a[0]);
   Assert.AreEqual(c[1],a[1]);
   Assert.AreEqual(c[2],a[2]);
   Assert.AreEqual(c[3],a[3]);
 }

 public void CtorCopy()
 {
   ComplexDoubleVector a = new ComplexDoubleVector(new double[2]{0,1});
   ComplexDoubleVector b = new ComplexDoubleVector(a);
   
   Assert.AreEqual(b.Length,a.Length);
   Assert.AreEqual(b[0],a[0]);
   Assert.AreEqual(b[1],a[1]);
 }

 public void Subtract()
 {
   ComplexDoubleVector a = new ComplexDoubleVector(new double[4]{0,1,2,3});
   ComplexDoubleVector b = new ComplexDoubleVector(new double[4]{4,5,6,7});
   ComplexDoubleVector c = new ComplexDoubleVector(a.Length);
   ComplexDoubleVector d = new ComplexDoubleVector(b.Length);
   
   c = a-b;
   d = ComplexDoubleVector.Subtract(a,b);
   
   Assert.AreEqual(c[0],a[0]-b[0]);
   Assert.AreEqual(c[1],a[1]-b[1]);
   Assert.AreEqual(c[2],a[2]-b[2]);
   Assert.AreEqual(c[3],a[3]-b[3]);
   
   Assert.AreEqual(d[0],c[0]);
   Assert.AreEqual(d[1],c[1]);
   Assert.AreEqual(d[2],c[2]);
   Assert.AreEqual(d[3],c[3]);
   
   a.Subtract(b);
   
   Assert.AreEqual(c[0],a[0]);
   Assert.AreEqual(c[1],a[1]);
   Assert.AreEqual(c[2],a[2]);
   Assert.AreEqual(c[3],a[3]);
 }

 public void Negate()
 {
   double[] vec = new double[4]{0,1,2,3};
   ComplexDoubleVector a = new ComplexDoubleVector(vec);
   ComplexDoubleVector b = -a;
   
   a = ComplexDoubleVector.Negate(a);
   
   Assert.AreEqual(-(Complex)vec[0],a[0]);
   Assert.AreEqual(-(Complex)vec[1],a[1]);
   Assert.AreEqual(-(Complex)vec[2],a[2]);
   Assert.AreEqual(-(Complex)vec[3],a[3]);
   
   Assert.AreEqual(-(Complex)vec[0],b[0]);
   Assert.AreEqual(-(Complex)vec[1],b[1]);
   Assert.AreEqual(-(Complex)vec[2],b[2]);
   Assert.AreEqual(-(Complex)vec[3],b[3]);
 }

 public void Axpy()
 {
   ComplexDoubleVector a = new ComplexDoubleVector(new double[4]{0,1,2,3});
   Double scal = 3;
   ComplexDoubleVector b = new ComplexDoubleVector(4);
   
   b.Axpy(scal,a);
   a.Scale(scal);
   
   Assert.AreEqual(a[0],b[0]);
   Assert.AreEqual(a[1],b[1]);
   Assert.AreEqual(a[2],b[2]);
   Assert.AreEqual(a[3],b[3]);
 }