public static float Compute(ComplexFloat a, ComplexFloat b)
 {
   ComplexFloat temp = a * ComplexMath.Conjugate(a);
   temp += (b * ComplexMath.Conjugate(b));
   float ret = ComplexMath.Absolute(temp);
   return (float)System.Math.Sqrt(ret);
 }

 internal static int Compute( int m, int n, int k, ComplexFloat[] A, int lda, ComplexFloat[] tau ){
   ArgumentCheck(m, n, k, A, lda, tau);
   if (tau.Length < System.Math.Max(1, k) ){
     throw new ArgumentException("tau must be at least max(1,k).");
   }
   
   return dna_lapack_cungqr(Configuration.BlockSize, m, n, k, A, lda, tau);
 }

 internal static int Compute( int m, int n, ComplexFloat[] A, int lda, out float[] d, out float[] e, out ComplexFloat[] tauq, out ComplexFloat[] taup ){
   ArgumentCheck(m, n, A, lda);
   d = new float[System.Math.Max(1, System.Math.Min(m,n))];
   e = new float[System.Math.Max(1, System.Math.Min(m,n))];
   tauq = new ComplexFloat[System.Math.Max(1, System.Math.Min(m,n))];
   taup = new ComplexFloat[System.Math.Max(1, System.Math.Min(m,n))];
   
   return dna_lapack_cgebrd(Configuration.BlockSize, m, n, A, lda, d, e, tauq, taup);
 }

 ///<summary>Constructor for <c>ComplexFloatVector</c> from <c>ComplexFloat</c> array</summary>
 ///<param name="values">Array of <c>ComplexFloat</c> to convert into <c>ComplexFloatVector</c>.</param>
 ///<exception cref="ArgumentNullException">Exception thrown if null passed as 'value' parameter.</exception>
 public ComplexFloatVector(ComplexFloat[] values)
 {
   if (values == null)
   {
     throw new ArgumentNullException("Array cannot be null");
   }
   data = new ComplexFloat[values.Length];
   for (int i = 0; i < values.Length; ++i)
   {
     data[i] = values[i];
   }
 }

 internal static int Compute( Vector vect, Side side, Transpose trans, int m, int n, int k, ComplexFloat[] A, int lda, ComplexFloat[] tau, ComplexFloat[] C, int ldc ){
   ArgumentCheck(vect,side, m, n, k, A, lda, tau, C, ldc);
   if( side == Side.Left){
     if (tau.Length < System.Math.Max(1, System.Math.Min(m,k))){
       throw new ArgumentException("tau must be at least max(1,k).");
     }
   }else{
     if (tau.Length < System.Math.Max(1, System.Math.Min(n,k))){
       throw new ArgumentException("tau must be at least max(1,k).");
     }
   }
   
   return dna_lapack_cunmbr(Configuration.BlockSize, vect, side, trans, m, n, k, A, lda, tau, C, ldc);
 }

    ///<summary>Compute the function of this class</summary>

    internal static ComplexFloat Compute( int n, ComplexFloat[] X, int incx, ComplexFloat[] Y, int incy )
    {
      if ( n < 0)
      { 
        return ComplexFloat.Zero;
      }
      ArgumentCheck(n, X, X.Length, ref incx, Y, Y.Length, ref incy);
      
      ComplexFloat ret = new ComplexFloat(0);
#if MANAGED
      int ix = 0;
      int iy = 0;
      for ( int i = 0; i < n; ++i ) 
      {
        ret += (ComplexMath.Conjugate(X[ix]) * Y[iy]);
        ix += incx;
        iy += incy;
      }
#else
      dna_blas_cdotc(n, X, incx, Y, incy, ref ret);
#endif
      return ret;
    }

    /// <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 ComplexFloatMatrix Solve(IROComplexFloatMatrix Y)
    {
      ComplexFloatMatrix 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
      ComplexFloat[] Inner;     // inner product
      ComplexFloat[] G;       // scaling constant
      ComplexFloat[] A;       // reference to current order coefficients
      ComplexFloat scalar;

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

      // setup zero order solution
      scalar = ComplexFloat.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++)
        {
//.........这里部分代码省略.........

    internal static int Compute( int n, ComplexFloat[] A, int lda, int[] ipiv ){
      ArgumentCheck(n,A,lda,ipiv);

      return dna_lapack_cgetri(Configuration.BlockSize, n,A,lda,ipiv);
    }

    /// <overloads>
    /// There are two permuations of the constructor, both require a parameter corresponding
    /// to the left-most column of a Toeplitz matrix.
    /// </overloads>
    /// <summary>
    /// Constructor with <c>ComplexFloatVector</c> parameter.
    /// </summary>
    /// <param name="T">The left-most column of the Toeplitz matrix.</param>
    /// <exception cref="ArgumentNullException">
    /// <B>T</B> is a null reference.
    /// </exception>
    /// <exception cref="RankException">
    /// The length of <B>T</B> is zero.
    /// </exception>
    public ComplexFloatSymmetricLevinson(IROComplexFloatVector T)
    {
      // check parameter
      if (T == null)
      {
        throw new System.ArgumentNullException("T");
      }
      else if (T.Length == 0)
      {
        throw new RankException("The length of T is zero.");
      }

      // save the vector
      m_LeftColumn = new ComplexFloatVector(T);
      m_Order = m_LeftColumn.Length;

      // allocate memory for lower triangular matrix
      m_LowerTriangle = new ComplexFloat[m_Order][];
      for (int i = 0; i < m_Order; i++)
      {
        m_LowerTriangle[i] = new ComplexFloat[i+1];
      }

      // allocate memory for diagonal
      m_Diagonal = new ComplexFloat[m_Order];

    }

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

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

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

      if (m_Order > 1)
      {
        // fill bottom row (reverse order)
        ComplexFloat[] 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)
      {
        ComplexFloat[] top = new ComplexFloat[m_Order];
        Array.Copy(m_LeftColumn.data, 0, top, 0, m_Order);
        tm.SetRow(0, top);

        // fill bottom row (reverse order)
        ComplexFloat[] bottom = new ComplexFloat[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++)
        {
          ComplexFloat[] temp = new ComplexFloat[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;
    }

 internal static int Compute( int m, int n, ComplexFloat[] a, float[] s, ComplexFloat[] u, ComplexFloat[] v  ){
   if( a == null ){
     throw new ArgumentNullException("a", "a cannot be null.");
   }
   return dna_lapack_cgesvd(m, n, a, m, s, u, m, v, n);
 }

		/// <summary>
		/// Invert a square Toeplitz matrix.
		/// </summary>
		/// <param name="col">The left-most column of the Toeplitz matrix.</param>
		/// <param name="row">The top-most row of the Toeplitz matrix.</param>
		/// <returns>The inverse matrix.</returns>
		/// <exception cref="ArgumentNullException">
		/// <B>col</B> is a null reference.
		/// <para>or</para>
		/// <para><B>row</B> is a null reference.</para>
		/// </exception>
		/// <exception cref="RankException">
		/// The length of <B>col</B> is 0,
		/// <para>or</para>
		/// <para>the lengths of <B>col</B> and <B>row</B> are not equal.</para>
		/// </exception>
		/// <exception cref="SingularMatrixException">
		/// The Toeplitz matrix or one of the the leading sub-matrices is singular.
		/// </exception>
		/// <exception cref="ArithmeticException">
		/// The values of the first element of <B>col</B> and <B>row</B> are not equal.
		/// </exception>
		/// <remarks>
		/// This static member combines the <b>UDL</b> decomposition and Trench's algorithm into a
		/// single algorithm. It requires minimal data storage, compared to the non-static member
		/// and suffers from no speed penalty in comparision.
		/// <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 ComplexFloatMatrix Inverse(IROComplexFloatVector col, IROComplexFloatVector row)
		{
			// check parameters
			if (col == null)
			{
				throw new System.ArgumentNullException("col");
			}
			else if (col.Length == 0)
			{
				throw new RankException("The length of col is zero.");
			}
			else if (row == null)
			{
				throw new System.ArgumentNullException("row");
			}
			else if (col.Length != row.Length)
			{
				throw new RankException("The lengths of col and row are not equal.");
			}
			else if (col[0] != row[0])
			{
				throw new ArithmeticException("The values of the first element of col and row are not equal.");
			}

			// check if leading diagonal is zero
			if (col[0] == ComplexFloat.Zero)
			{
				throw new SingularMatrixException("One of the leading sub-matrices is singular.");
			}

			// decompose matrix
			int order = col.Length;
			ComplexFloat[] A = new ComplexFloat[order];
			ComplexFloat[] B = new ComplexFloat[order];
			ComplexFloat[] Z = new ComplexFloat[order];
			ComplexFloat Q, S, Ke, Kr, e;
			int i, j, k, l;

			// setup the zero order solution
			A[0] = ComplexFloat.One;
			B[0] = ComplexFloat.One;
			e = ComplexFloat.One / col[0];

			for (i = 1; i < order; i++)
			{
				// calculate inner products
				Q = ComplexFloat.Zero;
				for (j = 0, l = 1; j < i; j++, l++)
				{
					Q += col[l] * A[j];
				}

				S = ComplexFloat.Zero;
				for (j = 0, l = 1; j < i; j++, l++)
				{
					S += row[l] * B[j];
				}

				// reflection coefficients
				Kr = -S * e;
				Ke = -Q * e;

				// update lower triangle (in temporary storage)
				Z[0] = ComplexFloat.Zero;
				Array.Copy(A, 0, Z, 1, i);
				for (j = 0, l = i - 1; j < i; j++, l--)
				{
					Z[j] += Ke * B[l];
				}
//.........这里部分代码省略.........

 public static int Compute( int n, int ilo, int ihi, ComplexFloat[] A, int lda, ComplexFloat[] tau ){
   ArgumentCheck(n, ilo, ihi, A, lda, tau);
   return dna_lapack_cgehrd(Configuration.BlockSize, n, ilo, ihi, A, lda, tau);
 }

 private static extern void dna_blas_ctrsm( Order order, Side side, UpLo uplo, Transpose transA, Diag diag, int m,int n, ref ComplexFloat alpha, [In]ComplexFloat[] A, int lda, [In,Out]ComplexFloat[] B, int ldb );

		internal static void Compute(int n, ComplexFloat[] X, int incx, ComplexFloat[] Y, int incy)
		{
			if (n < 0)
			{
				return;
			}
			ArgumentCheck(n, X, X.Length, ref incx, Y, Y.Length, ref incy);

#if MANAGED
			int ix = 0;
			int iy = 0;
			for (int i = 0; i < n; ++i)
			{
				Y[iy] = X[ix];
				ix += incx;
				iy += incy;
			}
#else
      dna_blas_ccopy(n, X, incx, Y, incy);
#endif
		}

 internal static int Compute( int m, int n, ComplexFloat[] A, int lda, int[] jpvt, out ComplexFloat[] tau  ){
   ArgumentCheck(m, n, A, lda, jpvt);
   tau = new ComplexFloat[System.Math.Max(1, System.Math.Min(m,n))];
   
   return dna_lapack_cgeqp3(m,n,A,lda,jpvt,tau);
 }

 private static extern void dna_blas_cdotc(int N, [In]ComplexFloat[] X, int incX, [In]ComplexFloat[] Y, int incY, ref ComplexFloat dotc);

		internal static void Compute(Order order, Transpose transA, int m, int n, ComplexFloat alpha, ComplexFloat[] A, int lda, ComplexFloat[] X, int incx, ComplexFloat beta, ComplexFloat[] Y, int incy)
		{
			ArgumentCheck(order, transA, m, n, A, A.Length, lda, X, X.Length, ref incx, Y, Y.Length, ref incy);

#if MANAGED
			int lenX = m;
			int lenY = n;

			if (transA == Transpose.NoTrans)
			{
				lenX = n;
				lenY = m;
			}
			if (beta == ComplexFloat.Zero)
			{
				for (int i = 0, iy = 0; i < lenY; ++i, iy += incy)
				{
					Y[iy] = 0;
				}
			}
			else if (beta != ComplexFloat.One)
			{
				for (int i = 0, iy = 0; i < lenY; ++i, iy += incy)
				{
					Y[iy] *= beta;
				}
			}
			if (alpha == ComplexFloat.Zero)
			{
				return;
			}

			ComplexFloat add_val = new ComplexFloat(0);
			if ((order == Order.RowMajor && transA == Transpose.NoTrans)
				|| (order == Order.ColumnMajor && transA == Transpose.Trans))
			{
				for (int i = 0, iy = 0; i < lenY; ++i, iy += incy)
				{
					add_val = 0;
					for (int j = 0, jx = 0; j < lenX; ++j, jx += incx)
					{
						add_val += X[jx] * A[i * lda + j];
					}
					Y[iy] += alpha * add_val;
				}
			}
			else if ((order == Order.RowMajor && transA == Transpose.Trans)
				|| (order == Order.ColumnMajor && transA == Transpose.NoTrans))
			{
				for (int j = 0, jx = 0; j < lenX; ++j, jx += incx)
				{
					add_val = alpha * X[jx];
					if (add_val != ComplexFloat.Zero)
					{
						for (int i = 0, iy = 0; i < lenY; ++i, iy += incy)
						{
							Y[iy] += add_val * A[j * lda + i];
						}
					}
				}
			}
			else if (order == Order.RowMajor && transA == Transpose.ConjTrans)
			{
				for (int j = 0, jx = 0; j < lenX; ++j, jx += incx)
				{
					add_val = alpha * X[jx];
					if (add_val != ComplexFloat.Zero)
					{
						for (int i = 0, iy = 0; i < lenY; ++i, iy += incy)
						{
							Y[iy] += add_val * ComplexMath.Conjugate(A[j * lda + i]);
						}
					}
				}
			}
			else
			{
				for (int i = 0, iy = 0; i < lenY; ++i, iy += incy)
				{
					add_val = 0;
					for (int j = 0, jx = 0; j < lenX; ++j, jx += incx)
					{
						add_val += X[jx] * ComplexMath.Conjugate(A[i * lda + j]);
					}
					Y[iy] += alpha * add_val;
				}
			}
#else
      dna_blas_cgemv(order, transA, m, n, ref alpha, A, lda, X, incx, ref beta, Y, incy);
#endif
		}

 private static extern void dna_blas_cgemv( Order order, Transpose TransA, int M, int N, ref ComplexFloat alpha, [In]ComplexFloat[] A, int lda, [In]ComplexFloat[] X, int incX, ref ComplexFloat beta, [In,Out]ComplexFloat[] Y, int incY);

    internal static float Compute(int n, ComplexFloat[] X, int incx)
    {
      if( n < 1 ) 
      {
        return 0;
      }
      if ( incx <= 0 ) 
      {
        return 0;
      }
      ArgumentCheck(n, X, incx, X.Length);
      
      float ret = 0;

#if MANAGED
      int ix = 0;
      double temp = 0;
      for ( int i = 0; i < n; ++i ) 
      {
        temp += System.Math.Abs(X[ix].Real) + System.Math.Abs(X[ix].Imag);
        ix += incx;
      }
      ret = (float)temp;
#else
      ret = dna_blas_scasum(n, X, incx);
#endif
      return ret;
    }