C# ComplexFloatMatrix类代码示例

		public void Current()
		{
			ComplexFloatMatrix test = new ComplexFloatMatrix(new ComplexFloat[2, 2] { { 1f, 2f }, { 3f, 4f } });
			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 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()
 {
   ComplexFloatMatrix test = new ComplexFloatMatrix(2,2);
   
   Assert.AreEqual(test.RowLength, 2);
   Assert.AreEqual(test.ColumnLength, 2);
   Assert.AreEqual(test[0,0], ComplexFloat.Zero);
   Assert.AreEqual(test[0,1], ComplexFloat.Zero);
   Assert.AreEqual(test[1,0], ComplexFloat.Zero);
   Assert.AreEqual(test[1,1], ComplexFloat.Zero);
 }

 public void CtorInitialValues()
 {
   ComplexFloatMatrix test = new ComplexFloatMatrix(2,2,new ComplexFloat(1,1));
   
   Assert.AreEqual(test.RowLength, 2);
   Assert.AreEqual(test.ColumnLength, 2);
   ComplexFloat value = new ComplexFloat(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);
 }

    ///<summary>Computes the algorithm.</summary>
    protected override void InternalCompute()
    {
#if MANAGED
      l = new ComplexFloatMatrix(matrix);
      for (int j = 0; j < order; j++) 
      {
        ComplexFloat[] rowj = l.data[j];
        float d = 0.0f;
        for (int k = 0; k < j; k++) 
        {
          ComplexFloat[] rowk = l.data[k];
          ComplexFloat s = ComplexFloat.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 ComplexFloat((float)System.Math.Sqrt(d));
        for (int k = j+1; k < order; k++) 
        {
          l.data[j][k] = ComplexFloat.Zero;
        }
      }
#else
            ComplexFloat[] factor = new ComplexFloat[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 ComplexFloatMatrix(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    
    }

		public static bool AreEqual(ComplexFloatMatrix f1, ComplexFloatMatrix 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 ComplexFloatLUDecompTest() 
 {
   ComplexFloatMatrix a = new ComplexFloatMatrix(3);
   a[0,0] = new ComplexFloat(-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 ComplexFloatLUDecomp(a);
 }

 static ComplexFloatCholeskyDecompTest() 
 {
   ComplexFloatMatrix a = new ComplexFloatMatrix(3);
   a[0,0] = 2;
   a[0,1] = new ComplexFloat(1,-1);
   a[0,2] = 0;
   a[1,0] = new ComplexFloat(1,-1);
   a[1,1] = 2;
   a[1,2] = 0;
   a[2,0] = 0;
   a[2,1] = 0;
   a[2,2] = 3;
   cd = new ComplexFloatCholeskyDecomp(a);
 }

    public void SquareDecomp()
    {
      ComplexFloatMatrix a = new ComplexFloatMatrix(3);
      a[0,0] = new ComplexFloat(1.1f, 1.1f);
      a[0,1] = new ComplexFloat(2.2f, -2.2f);
      a[0,2] = new ComplexFloat(3.3f, 3.3f);
      a[1,0] = new ComplexFloat(4.4f, -4.4f);
      a[1,1] = new ComplexFloat(5.5f, 5.5f);
      a[1,2] = new ComplexFloat(6.6f, -6.6f);
      a[2,0] = new ComplexFloat(7.7f, 7.7f);
      a[2,1] = new ComplexFloat(8.8f, -8.8f);
      a[2,2] = new ComplexFloat(9.9f, 9.9f);
      
      ComplexFloatQRDecomp qrd = new ComplexFloatQRDecomp(a);
      ComplexFloatMatrix qq = qrd.Q.GetConjugateTranspose()*qrd.Q;
      ComplexFloatMatrix qr = qrd.Q*qrd.R;
      ComplexFloatMatrix I = ComplexFloatMatrix.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-6);
      
      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 < 2.4E-6);
    }

    ///<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 ComplexFloatCholeskyDecomp(IROComplexFloatMatrix 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 ComplexFloatMatrix(matrix);
    }

 public void CtorCopy()
 {
   ComplexFloatMatrix a = new ComplexFloatMatrix(2,2);
   a[0,0] = new ComplexFloat(1,1);     
   a[0,1] = new ComplexFloat(2,2);     
   a[1,0] = new ComplexFloat(3,3);     
   a[1,1] = new ComplexFloat(4,4); 
   
   ComplexFloatMatrix b = new ComplexFloatMatrix(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()
 {
   ComplexFloatMatrix b = new ComplexFloatMatrix(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;
   ComplexFloatCholeskyDecomp dcd = new ComplexFloatCholeskyDecomp(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 ComplexFloatMatrix(3);
   a[0,0] = new ComplexFloat(1.1f, 1.1f);
   a[0,1] = new ComplexFloat(2.2f, -2.2f);
   a[0,2] = new ComplexFloat(3.3f, 3.3f);
   a[1,0] = new ComplexFloat(4.4f, -4.4f);
   a[1,1] = new ComplexFloat(5.5f, 5.5f);
   a[1,2] = new ComplexFloat(6.6f, -6.6f);
   a[2,0] = new ComplexFloat(7.7f, 7.7f);
   a[2,1] = new ComplexFloat(8.8f, -8.8f);
   a[2,2] = new ComplexFloat(9.9f, 9.9f);
   svd = new ComplexFloatSVDDecomp(a, true);
   
   wa = new ComplexFloatMatrix(2,4);
   wa[0,0] = new ComplexFloat(1.1f, 1.1f);
   wa[0,1] = new ComplexFloat(2.2f, -2.2f);
   wa[0,2] = new ComplexFloat(3.3f, 3.3f);
   wa[0,3] = new ComplexFloat(4.4f, -4.4f);
   wa[1,0] = new ComplexFloat(5.5f, 5.5f);
   wa[1,1] = new ComplexFloat(6.6f, -6.6f);
   wa[1,2] = new ComplexFloat(7.7f, 7.7f);
   wa[1,3] = new ComplexFloat(8.8f, -8.8f);
   wsvd = new ComplexFloatSVDDecomp(wa, true);
     
   la = new ComplexFloatMatrix(4,2);
   la[0,0] = new ComplexFloat(1.1f, 1.1f);
   la[0,1] = new ComplexFloat(2.2f, -2.2f);
   la[1,0] = new ComplexFloat(3.3f, 3.3f);
   la[1,1] = new ComplexFloat(4.4f, -4.4f);
   la[2,0] = new ComplexFloat(5.5f, 5.5f);
   la[2,1] = new ComplexFloat(6.6f, -6.6f);
   la[3,0] = new ComplexFloat(7.7f, 7.7f);
   la[3,1] = new ComplexFloat(8.8f, -8.8f);
   lsvd = new ComplexFloatSVDDecomp(la, true);
 } 

    /// <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 ComplexFloatMatrix Inverse(IROComplexFloatVector T)
    {

      ComplexFloatMatrix 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 ComplexFloatMatrix(1);
        X[0, 0] = ComplexFloat.One / T[0];
      }
      else
      {

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

        // calculate the predictor coefficients
        ComplexFloatVector Y = ComplexFloatSymmetricLevinson.YuleWalker(T);

        // calculate gamma
        f = T[0];
        for (i = 1, j = 0; i < N; i++, j++)
        {
          f += T[i] * Y[j];
        }
        g = ComplexFloat.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>
    /// Solve a symmetric square Toeplitz system with a right-side matrix.
    /// </summary>
    /// <param name="T">The left-most column of the Toeplitz matrix.</param>
    /// <param name="Y">The right-side matrix of the system.</param>
    /// <returns>The solution matrix.</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 number of rows in <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
    /// matrix and <B>Y</B> is a known matrix.
    /// <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 ComplexFloatMatrix Solve(IROComplexFloatVector T, IROComplexFloatMatrix Y)
    {

      ComplexFloatMatrix 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.Columns)
      {
        throw new RankException("The length of T and Y are not equal.");
      }
      else
      {

        // allocate memory
        int N = T.Length;
        int M = Y.Rows;
        X = new ComplexFloatMatrix(N, M);                 // solution matrix
        ComplexFloatVector Z = new ComplexFloatVector(N);       // temporary storage vector
        ComplexFloat e;                                   // prediction error
        int i, j, l, m;

        // 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.");
        }
        for (m = 0; m < M; m++)
        {
          X[0, m] = Y[0,m] / T[0];
        }

        if (N > 1)
        {

          ComplexFloatVector a = new ComplexFloatVector(N - 1);   // prediction coefficients
          ComplexFloat p;                                   // reflection coefficient
          ComplexFloat inner;                               // inner product
          ComplexFloat k;

          // 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
            p = -(inner / e);
            for (j = 0, l = i - 2; j < i - 1; j++, l--)
            {
              Z[j] = a[j] + p * a[l];
            }

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

            a[i - 1] = p;
            e *= (ComplexFloat.One - p * p);

//.........这里部分代码省略.........