public void Current()
		{
			ComplexDoubleMatrix test = new ComplexDoubleMatrix(new Complex[2, 2] { { 1, 2 }, { 3, 4 } });
			IEnumerator enumerator = test.GetEnumerator();
			bool movenextresult;

			movenextresult = enumerator.MoveNext();
			Assert.IsTrue(movenextresult);
			Assert.AreEqual(enumerator.Current, test[0, 0]);

			movenextresult = enumerator.MoveNext();
			Assert.IsTrue(movenextresult);
			Assert.AreEqual(enumerator.Current, test[1, 0]);

			movenextresult = enumerator.MoveNext();
			Assert.IsTrue(movenextresult);
			Assert.AreEqual(enumerator.Current, test[0, 1]);

			movenextresult = enumerator.MoveNext();
			Assert.IsTrue(movenextresult);
			Assert.AreEqual(enumerator.Current, test[1, 1]);

			movenextresult = enumerator.MoveNext();
			Assert.IsFalse(movenextresult);
		}

    /// <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 CtorDimensions()
		{
			ComplexDoubleMatrix test = new ComplexDoubleMatrix(2, 2);

			Assert.AreEqual(test.RowLength, 2);
			Assert.AreEqual(test.ColumnLength, 2);
			Assert.AreEqual(test[0, 0], Complex.Zero);
			Assert.AreEqual(test[0, 1], Complex.Zero);
			Assert.AreEqual(test[1, 0], Complex.Zero);
			Assert.AreEqual(test[1, 1], Complex.Zero);
		}

		public void CtorInitialValues()
		{
			ComplexDoubleMatrix test = new ComplexDoubleMatrix(2, 2, new Complex(1, 1));

			Assert.AreEqual(test.RowLength, 2);
			Assert.AreEqual(test.ColumnLength, 2);
			Complex value = new Complex(1, 1);
			Assert.AreEqual(test[0, 0], value);
			Assert.AreEqual(test[0, 1], value);
			Assert.AreEqual(test[1, 0], value);
			Assert.AreEqual(test[1, 1], value);
		}

 public static bool AreEqual(ComplexDoubleMatrix f1, ComplexDoubleMatrix f2, float delta) 
 {
   if (f1.RowLength != f2.RowLength) return false;
   if (f1.ColumnLength != f2.ColumnLength) return false;
   for (int i = 0; i < f1.RowLength; i++) 
   {
     for (int j = 0; j < f1.ColumnLength; j++) 
     {
       if (!AreEqual(f1[i, j], f2[i, j], delta))
         return false;
     }
   }
   return true;
 }

 static ComplexDoubleCholeskyDecompTest()
 {
   ComplexDoubleMatrix a = new ComplexDoubleMatrix(3);
   a[0,0] = 2;
   a[0,1] = new Complex(1,-1);
   a[0,2] = 0;
   a[1,0] = new Complex(1,-1);
   a[1,1] = 2;
   a[1,2] = 0;
   a[2,0] = 0;
   a[2,1] = 0;
   a[2,2] = 3;
   cd = new ComplexDoubleCholeskyDecomp(a);
 }

 static ComplexDoubleLUDecompTest() 
 {
   ComplexDoubleMatrix a = new ComplexDoubleMatrix(3);
   a[0,0] = new Complex(-1,1);
   a[0,1] = 5;
   a[0,2] = 6;
   a[1,0] = 3;
   a[1,1] = -6;
   a[1,2] = 1;
   a[2,0] = 6;
   a[2,1] = 8;
   a[2,2] = 9;
   lu = new ComplexDoubleLUDecomp(a);
 }

    ///<summary>Computes the algorithm.</summary>
    protected override void InternalCompute()
    {  
#if MANAGED
      l = new ComplexDoubleMatrix(matrix);
      for (int j = 0; j < order; j++) 
      {
        Complex[] rowj = l.data[j];
        double d = 0.0;
        for (int k = 0; k < j; k++) 
        {
          Complex[] rowk = l.data[k];
          Complex s = Complex.Zero;
          for (int i = 0; i < k; i++) 
          {
            s += rowk[i]*rowj[i];
          }
          rowj[k] = s = (matrix.data[j][k] - s)/l.data[k][k];
          d = d + (s*ComplexMath.Conjugate(s)).Real;
        }
        d = matrix.data[j][j].Real - d;
        if ( d <= 0.0 ) 
        {
          ispd = false;
          return;
        }
        l.data[j][j] = new Complex(System.Math.Sqrt(d));
        for (int k = j+1; k < order; k++) 
        {
          l.data[j][k] = Complex.Zero;
        }
      }
#else
            Complex[] factor = new Complex[matrix.data.Length];
            Array.Copy(matrix.data, factor, matrix.data.Length);
            int status = Lapack.Potrf.Compute(Lapack.UpLo.Lower, order, factor, order);
            if (status != 0 ) {
                ispd = false;
            }
            l = new ComplexDoubleMatrix(order);
            l.data = factor;
            for (int i = 0; i < order; i++) {
                for (int j = 0; j < order; j++) {
                    if ( j > i) {
                        l.data[j*order+i] = 0;
                    }
                }
            }
#endif    
    }

    ///<summary>Constructor for Cholesky decomposition class. The constructor performs the factorization of a Hermitian positive
    ///definite matrax and the Cholesky factored matrix is accessible by the <c>Factor</c> property. The factor is the lower 
    ///triangular factor.</summary>
    ///<param name="matrix">The matrix to factor.</param>
    ///<exception cref="ArgumentNullException">matrix is null.</exception>
    ///<exception cref="NotSquareMatrixException">matrix is not square.</exception>
    ///<remarks>This class only uses the lower triangle of the input matrix. It ignores the
    ///upper triangle.</remarks>
    public ComplexDoubleCholeskyDecomp(IROComplexDoubleMatrix matrix)
    {
      if ( matrix == null ) 
      {
        throw new System.ArgumentNullException("matrix cannot be null.");
      }

      if ( matrix.Rows != matrix.Columns ) 
      {
        throw new NotSquareMatrixException("Matrix must be square.");
      }

      order = matrix.Columns;
      this.matrix = new ComplexDoubleMatrix(matrix);
    }

		public void SquareDecomp()
		{
			ComplexDoubleMatrix a = new ComplexDoubleMatrix(3);
			a[0, 0] = new Complex(1.1, 1.1);
			a[0, 1] = new Complex(2.2, -2.2);
			a[0, 2] = new Complex(3.3, 3.3);
			a[1, 0] = new Complex(4.4, -4.4);
			a[1, 1] = new Complex(5.5, 5.5);
			a[1, 2] = new Complex(6.6, -6.6);
			a[2, 0] = new Complex(7.7, 7.7);
			a[2, 1] = new Complex(8.8, -8.8);
			a[2, 2] = new Complex(9.9, 9.9);

			ComplexDoubleQRDecomp qrd = new ComplexDoubleQRDecomp(a);
			ComplexDoubleMatrix qq = qrd.Q.GetConjugateTranspose() * qrd.Q;
			ComplexDoubleMatrix qr = qrd.Q * qrd.R;
			ComplexDoubleMatrix I = ComplexDoubleMatrix.CreateIdentity(3);

			// determine the maximum relative error
			double MaxError = 0.0;
			for (int i = 0; i < 3; i++)
			{
				for (int j = 0; i < 3; i++)
				{
					double E = ComplexMath.Absolute((qq[i, j] - I[i, j]));
					if (E > MaxError)
					{
						MaxError = E;
					}
				}
			}
			Assert.IsTrue(MaxError < 1.0E-14);

			MaxError = 0.0;
			for (int i = 0; i < 3; i++)
			{
				for (int j = 0; i < 3; i++)
				{
					double E = ComplexMath.Absolute((qr[i, j] - a[i, j]) / a[i, j]);
					if (E > MaxError)
					{
						MaxError = E;
					}
				}
			}
			Assert.IsTrue(MaxError < 1.0E-14);
		}

		public void CtorCopy()
		{
			ComplexDoubleMatrix a = new ComplexDoubleMatrix(2, 2);
			a[0, 0] = new Complex(1, 1);
			a[0, 1] = new Complex(2, 2);
			a[1, 0] = new Complex(3, 3);
			a[1, 1] = new Complex(4, 4);

			ComplexDoubleMatrix b = new ComplexDoubleMatrix(a);

			Assert.AreEqual(a.RowLength, b.RowLength);
			Assert.AreEqual(a.ColumnLength, b.ColumnLength);
			Assert.AreEqual(a[0, 0], b[0, 0]);
			Assert.AreEqual(a[0, 1], b[0, 1]);
			Assert.AreEqual(a[1, 0], b[1, 0]);
			Assert.AreEqual(a[1, 1], b[1, 1]);
		}

 public void NonSymmFactorTest()
 {
   ComplexDoubleMatrix b = new ComplexDoubleMatrix(3);
   b[0,0] = 2;
   b[0,1] = 1;
   b[0,2] = 1;
   b[1,0] = 1;
   b[1,1] = 2;
   b[1,2] = 0;
   b[2,0] = 0;
   b[2,1] = 0;
   b[2,2] = 3;
   ComplexDoubleCholeskyDecomp dcd = new ComplexDoubleCholeskyDecomp(b);
   Assert.AreEqual(dcd.Factor[0,0].Real,1.414,TOLERENCE);
   Assert.AreEqual(dcd.Factor[0,1].Real,0.000,TOLERENCE);
   Assert.AreEqual(dcd.Factor[0,2].Real,0.000,TOLERENCE);
   Assert.AreEqual(dcd.Factor[1,0].Real,0.707,TOLERENCE);
   Assert.AreEqual(dcd.Factor[1,1].Real,1.225,TOLERENCE);
   Assert.AreEqual(dcd.Factor[1,2].Real,0.000,TOLERENCE);
   Assert.AreEqual(dcd.Factor[2,0].Real,0.000,TOLERENCE);
   Assert.AreEqual(dcd.Factor[2,1].Real,0.000,TOLERENCE);
   Assert.AreEqual(dcd.Factor[2,2].Real,1.732,TOLERENCE);
 }

		public void SetupTestCases()
		{
			a = new ComplexDoubleMatrix(3);
			a[0, 0] = new Complex(1.1, 1.1);
			a[0, 1] = new Complex(2.2, -2.2);
			a[0, 2] = new Complex(3.3, 3.3);
			a[1, 0] = new Complex(4.4, -4.4);
			a[1, 1] = new Complex(5.5, 5.5);
			a[1, 2] = new Complex(6.6, -6.6);
			a[2, 0] = new Complex(7.7, 7.7);
			a[2, 1] = new Complex(8.8, -8.8);
			a[2, 2] = new Complex(9.9, 9.9);
			svd = new ComplexDoubleSVDDecomp(a, true);

			wa = new ComplexDoubleMatrix(2, 4);
			wa[0, 0] = new Complex(1.1, 1.1);
			wa[0, 1] = new Complex(2.2, -2.2);
			wa[0, 2] = new Complex(3.3, 3.3);
			wa[0, 3] = new Complex(4.4, -4.4);
			wa[1, 0] = new Complex(5.5, 5.5);
			wa[1, 1] = new Complex(6.6, -6.6);
			wa[1, 2] = new Complex(7.7, 7.7);
			wa[1, 3] = new Complex(8.8, -8.8);
			wsvd = new ComplexDoubleSVDDecomp(wa, true);

			la = new ComplexDoubleMatrix(4, 2);
			la[0, 0] = new Complex(1.1, 1.1);
			la[0, 1] = new Complex(2.2, -2.2);
			la[1, 0] = new Complex(3.3, 3.3);
			la[1, 1] = new Complex(4.4, -4.4);
			la[2, 0] = new Complex(5.5, 5.5);
			la[2, 1] = new Complex(6.6, -6.6);
			la[3, 0] = new Complex(7.7, 7.7);
			la[3, 1] = new Complex(8.8, -8.8);
			lsvd = new ComplexDoubleSVDDecomp(la, true);
		}

    /// <summary>
    /// Solve a symmetric square Toeplitz system with a right-side matrix.
    /// </summary>
    /// <param name="Y">The right-hand side of the system.</param>
    /// <returns>The solution matrix.</returns>
    /// <exception cref="ArgumentNullException">
    /// Parameter <B>Y</B> is a null reference.
    /// </exception>
    /// <exception cref="RankException">
    /// The number of rows in <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
    /// a symmetric square Toeplitz matrix, <B>X</B> is the unknown solution matrix
    /// and <B>Y</B> is a known matrix.
    /// <para>
    /// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
    /// using the Levinson algorithm, and then calculates the solution matrix.
    /// </para>
    /// </remarks>
    public ComplexDoubleMatrix Solve(IROComplexDoubleMatrix Y)
    {
      ComplexDoubleMatrix X;

      // check parameters
      if (Y == null)
      {
        throw new System.ArgumentNullException("Y");
      }
      else if (m_Order != Y.Columns)
      {
        throw new RankException("The numer of rows in 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 M = Y.Rows;
      int i, j, l, m;     // index/loop variables
      Complex[] Inner;      // inner product
      Complex[] G;        // scaling constant
      Complex[] A;        // reference to current order coefficients
      Complex scalar;

      // allocate memory for solution
      X = new ComplexDoubleMatrix(m_Order, M);
      Inner = new Complex[M];
      G = new Complex[M];

      // setup zero order solution
      scalar = Complex.One / m_LeftColumn[0];
      for (m = 0; m < M; m++)
      {
#if MANAGED
        X.data[0][m] = scalar * Y[0,m];
#else

        X.data[m*m_Order] = scalar * Y[0,m];
#endif
      }

      // solve systems of increasing order
      for (i = 1; i < m_Order; i++)
      {
        // calculate inner product
        for (m = 0; m < M; m++)
        {
#if MANAGED
          Inner[m] = Y[i,m];
#else
          Inner[m] = Y[i,m];
#endif
        }

        for (j = 0, l = i; j < i; j++, l--)
        {
          scalar = m_LeftColumn[l];
          for (m = 0; m < M; m++)
          {
#if MANAGED
            Inner[m] -= scalar * X.data[j][m];
#else
            Inner[m] -= scalar * X.data[m*m_Order+j];
#endif
          }
        }

        // get the current predictor coefficients row
        A = m_LowerTriangle[i];

        // update the solution matrix
        for (m = 0; m < M; m++)
        {
//.........这里部分代码省略.........

    /// <summary>
    /// Get a copy of the Toeplitz matrix.
    /// </summary>
    public ComplexDoubleMatrix GetMatrix()
    {
      int i, j;

      // allocate memory for the matrix
      ComplexDoubleMatrix tm = new ComplexDoubleMatrix(m_Order);

#if MANAGED
      // fill top row
      Complex[] top = tm.data[0];
      Array.Copy(m_LeftColumn.data, 0, top, 0, m_Order);

      if (m_Order > 1)
      {
        // fill bottom row (reverse order)
        Complex[] bottom = tm.data[m_Order - 1];

        for (i = 0, j = m_Order - 1; i < m_Order; i++, j--)
        {
          bottom[i] = m_LeftColumn[j];
        }

        // fill rows in-between
        for (i = 1, j = m_Order - 1 ; j > 1; i++)
        {
          Array.Copy(top, 0, tm.data[i], i, j--);
          Array.Copy(bottom, j, tm.data[i], 0, i);
        }
      }
#else
      if (m_Order > 1)
      {
        Complex[] top = new Complex[m_Order];
        Array.Copy(m_LeftColumn.data, 0, top, 0, m_Order);
        tm.SetRow(0, top);

        // fill bottom row (reverse order)
        Complex[] bottom = new Complex[m_Order];

        for (i = 0, j = m_Order - 1; i < m_Order; i++, j--)
        {
          bottom[i] = m_LeftColumn[j];
        }

        // fill rows in-between
        for (i = 1, j = m_Order - 1 ; j > 0; i++)
        {
          Complex[] temp = new Complex[m_Order];
          Array.Copy(top, 0, temp, i, j--);
          Array.Copy(bottom, j, temp, 0, i);
          tm.SetRow(i, temp);
        }
      }
      else
      {
        Array.Copy(m_LeftColumn.data, 0, tm.data, 0, m_Order);
      }
#endif

      return tm;
    }

    /// <summary>
    /// Invert a symmetric square Toeplitz matrix.
    /// </summary>
    /// <param name="T">The left-most column of the symmetric Toeplitz matrix.</param>
    /// <returns>The inverse matrix.</returns>
    /// <exception cref="ArgumentNullException">
    /// <B>T</B> is a null reference.
    /// </exception>
    /// <exception cref="RankException">
    /// The length of <B>T</B> must be greater than zero.
    /// </exception>
    /// <exception cref="SingularMatrixException">
    /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
    /// </exception>
    /// <remarks>
    /// This static member combines the <b>UDL</b> decomposition and Trench's algorithm into a
    /// single algorithm. When compared to the non-static member it requires minimal data storage
    /// and suffers from no speed penalty.
    /// <para>
    /// Trench's algorithm requires <b>N</b> squared FLOPS, compared to <b>N</b> cubed FLOPS
    /// if we simply solved a linear Toeplitz system with a right-side identity matrix (<b>N</b> is the matrix order).
    /// </para>
    /// </remarks>
    public static ComplexDoubleMatrix Inverse(IROComplexDoubleVector T)
    {

      ComplexDoubleMatrix X;

      // check parameters
      if (T == null)
      {
        throw new System.ArgumentNullException("T");
      }
      else if (T.Length < 1)
      {
        throw new System.RankException("The length of T must be greater than zero.");
      }
      else if (T.Length == 1)
      {
        X = new ComplexDoubleMatrix(1);
        X[0, 0] = Complex.One / T[0];
      }
      else
      {

        int N = T.Length;
        Complex f, g;
        int i, j, l, k, m, n;
        X = new ComplexDoubleMatrix(N);

        // calculate the predictor coefficients
        ComplexDoubleVector Y = ComplexDoubleSymmetricLevinson.YuleWalker(T);

        // calculate gamma
        f = T[0];
        for (i = 1, j = 0; i < N; i++, j++)
        {
          f += T[i] * Y[j];
        }
        g = Complex.One / f;

        // calculate first row of inverse
        X[0, 0] = g;
        for (i = 1, j = 0; i < N; i++, j++)
        {
          X[0, i] = g * Y[j];
        }

        // calculate successive rows of upper wedge
        for (i = 0, j = 1, k = N - 2; i < N / 2; i++, j++, k--)
        {
          for (l = j, m = i, n = N-1-j; l < N - j; l++, m++, n--)
          {
            X[j, l] = X[i, m] + g * (Y[i] * Y[m] - Y[k] * Y[n]);
          }
        }

        // this is symmetric matrix ...
        for (i = 0; i <= N / 2; i++)
        {
          for (j = i + 1; j < N - i; j++)
          {
            X[j, i] = X[i, j];
          }
        }

        // and a persymmetric matrix.
        for (i = 0, j = N - 1; i < N; i++, j--)
        {
          for (k = 0, l = N - 1; k < j; k++, l--)
          {
            X[l, j] = X[i, k];
          }
        }

      }

      return X;
    }

 ///<summary> Constructor </summary>
 public ComplexDoubleMatrixEnumerator (ComplexDoubleMatrix matrix) 
 {
   m=matrix;
   index=-1;
   length=m.RowLength*m.ColumnLength;
 }

 public void ExplicitToComplexDoubleMatrix()
 {
   ComplexDoubleMatrix a = new ComplexDoubleMatrix(2,2);
   a[0,0] = 1;
   a[0,1] = 2;
   a[1,0] = 3;
   a[1,1] = 4;
   
   ComplexFloatMatrix b = ComplexFloatMatrix.ToComplexFloatMatrix(a);
   Assert.AreEqual(a.RowLength, b.RowLength);
   Assert.AreEqual(a.ColumnLength, b.ColumnLength);
   Assert.AreEqual(b[0,0], a[0,0]);
   Assert.AreEqual(b[0,1], a[0,1]);
   Assert.AreEqual(b[1,0], a[1,0]);
   Assert.AreEqual(b[1,1], a[1,1]);
 }

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

    private static void zdrot (ComplexDoubleMatrix A, int Col1, int Col2, double c, double s)
    {
      // applies a plane rotation, where the c=cos and s=sin
      // Col1, Col2 - rotated columns of A

      Complex z;
      for(int i=0; i < A.RowLength; i++)
      {
        z=c*A[i,Col1]+s*A[i,Col2];
        A[i,Col2]=c*A[i,Col2]-s*A[i,Col1];
        A[i,Col1]=z;
      }
    }