C# ComplexDoubleMatrix类代码示例

		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;
    }