C# FloatMatrix类代码示例

 static FloatQRDecompTest()
 {   
   FloatMatrix a = new FloatMatrix(3);
   a[0,0] = -1.0f;
   a[0,1] = 5.0f;
   a[0,2] = 6.0f;
   a[1,0] = 3.0f;
   a[1,1] = -6.0f;
   a[1,2] = 1.0f;
   a[2,0] = 6.0f;
   a[2,1] = 8.0f;
   a[2,2] = 9.0f;
   qr = new FloatQRDecomp(a);
   
   a = new FloatMatrix(2,3);
   a[0,0] = -1.0f;
   a[0,1] = 5.0f;
   a[0,2] = 6.0f;
   a[1,0] = 3.0f;
   a[1,1] = -6.0f;
   a[1,2] = 1.0f;
   wqr = new FloatQRDecomp(a);
   
   a = new FloatMatrix(3,2);
   a[0,0] = -1.0f;
   a[0,1] = 5.0f;
   a[1,0] = 3.0f;
   a[1,1] = -6.0f;
   a[2,0] = 6.0f;
   a[2,1] = 8.0f;
   lqr = new FloatQRDecomp(a);
 }

		public void Current()
		{
			FloatMatrix test = new FloatMatrix(new float[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.RowLength;
      int n = matrix.ColumnLength;

#if MANAGED
      int minmn = m < n ? m : n;
      r_ = new FloatMatrix(matrix); // create a copy
      FloatVector[] u = new FloatVector[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_ = FloatMatrix.CreateIdentity(m);
      for (int i = minmn - 1; i >= 0; i--) 
      {
        Householder.UA(u[i], q_, i, m - 1, i, m - 1);
      }
#else
      qr = new float[matrix.data.Length];
      Array.Copy(matrix.data, qr, matrix.data.Length);
      jpvt = new int[n];
      jpvt[0] = 1;
      Lapack.Geqp3.Compute(m, n, qr, m, jpvt, out tau);
      r_ = new FloatMatrix(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] = 0.0f;
          }
        }
      }
      q_ = new FloatMatrix(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] = 0.0f;
        }
      }

      if( m < n ){
        Lapack.Orgqr.Compute(m, m, m, q_.data, m, tau);
      } else{
        Lapack.Orgqr.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()
		{
			FloatMatrix test = new FloatMatrix(2, 2);
			Assert.AreEqual(test.RowLength, 2);
			Assert.AreEqual(test.ColumnLength, 2);
			Assert.AreEqual(test[0, 0], 0);
			Assert.AreEqual(test[0, 1], 0);
			Assert.AreEqual(test[1, 0], 0);
			Assert.AreEqual(test[1, 1], 0);
		}

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

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

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

#else
            float[] factor = new float[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 FloatMatrix(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
		}

 static FloatCholeskyDecompTest() 
 {
   FloatMatrix a = new FloatMatrix(3);
   a[0,0] = 2;
   a[0,1] = 1;
   a[0,2] = 0;
   a[1,0] = 1;
   a[1,1] = 2;
   a[1,2] = 0;
   a[2,0] = 0;
   a[2,1] = 0;
   a[2,2] = 3;
   cd = new FloatCholeskyDecomp(a);
 }

		static FloatLUDecompTest()
		{
			FloatMatrix a = new FloatMatrix(3);
			a[0, 0] = -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 FloatLUDecomp(a);
		}

		///<summary>Constructor for Cholesky decomposition class. The constructor performs the factorization of a symmetric 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 FloatCholeskyDecomp(IROFloatMatrix 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 FloatMatrix(matrix);
		}

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

			FloatMatrix b = new FloatMatrix(a);

			Assert.AreEqual(a.RowLength, b.RowLength);
			Assert.AreEqual(a.ColumnLength, b.ColumnLength);
			Assert.AreEqual(a[0, 0], a[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()
 {
   FloatMatrix b = new FloatMatrix(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;
   FloatCholeskyDecomp dcd = new FloatCholeskyDecomp(b);
   Assert.AreEqual(dcd.Factor[0,0],1.414,TOLERENCE);
   Assert.AreEqual(dcd.Factor[0,1],0.000,TOLERENCE);
   Assert.AreEqual(dcd.Factor[0,2],0.000,TOLERENCE);
   Assert.AreEqual(dcd.Factor[1,0],0.707,TOLERENCE);
   Assert.AreEqual(dcd.Factor[1,1],1.225,TOLERENCE);
   Assert.AreEqual(dcd.Factor[1,2],0.000,TOLERENCE);
   Assert.AreEqual(dcd.Factor[2,0],0.000,TOLERENCE);
   Assert.AreEqual(dcd.Factor[2,1],0.000,TOLERENCE);
   Assert.AreEqual(dcd.Factor[2,2],1.732,TOLERENCE);
 }

 public void SetupTestCases() 
 {
   a = new FloatMatrix(3);
   a[0,0] = 1.91f;
   a[0,1] = 9.82f;
   a[0,2] = 2.73f;
   a[1,0] = 8.64f;
   a[1,1] = 3.55f;
   a[1,2] = 7.46f;
   a[2,0] = 4.37f;
   a[2,1] = 6.28f;
   a[2,2] = 5.19f;
   svd = new FloatSVDDecomp(a, true);
   
   wa = new FloatMatrix(2,4);
   wa[0,0] = 1.91f;
   wa[0,1] = 9.82f;
   wa[0,2] = 2.73f;
   wa[0,3] = 8.64f;
   wa[1,0] = 3.55f;
   wa[1,1] = 7.46f;
   wa[1,2] = 4.37f;
   wa[1,3] = 6.28f;
   wsvd = new FloatSVDDecomp(wa, true);
     
   la = new FloatMatrix(4,2);
   la[0,0] = 1.91f;
   la[0,1] = 9.82f;
   la[1,0] = 2.73f;
   la[1,1] = 8.64f;
   la[2,0] = 3.55f;
   la[2,1] = 7.46f;
   la[3,0] = 4.37f;
   la[3,1] = 6.28f;
   lsvd = new FloatSVDDecomp(la, true);
 } 

    /// <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 FloatMatrix Solve(IROFloatVector T, IROFloatMatrix Y)
    {

      FloatMatrix 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 FloatMatrix(N, M);                 // solution matrix
        FloatVector Z = new FloatVector(N);       // temporary storage vector
        float e;                                   // prediction error
        int i, j, l, m;

        // setup zero order solution
        e = T[0];
        if (e == 0.0f)
        {
          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)
        {

          FloatVector a = new FloatVector(N - 1);   // prediction coefficients
          float p;                                   // reflection coefficient
          float inner;                               // inner product
          float 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 *= (1.0f - p * p);

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

    /// <summary>Finds the least squares solution of <c>A*X = B</c>, where <c>m &gt;= n</c></summary>
    /// <param name="B">A matrix with as many rows as A and any number of columns.</param>
    /// <returns>X that minimizes the two norm of <c>Q*R*X-B</c>.</returns>
    /// <exception cref="ArgumentException">Matrix row dimensions must agree.</exception>
    /// <exception cref="InvalidOperationException">Matrix is rank deficient or <c>m &lt; n</c>.</exception>
    public FloatMatrix Solve (IROFloatMatrix B) 
    {
      if (B.Rows != matrix.RowLength) 
      {
        throw new ArgumentException("Matrix row dimensions must agree.");
      }
      if (matrix.RowLength < matrix.ColumnLength) 
      {
        throw new System.InvalidOperationException("A must have at lest as a many rows as columns.");
      } 
      Compute();
      if (!this.isFullRank) 
      {
        throw new System.InvalidOperationException("Matrix is rank deficient.");
      }
      
      // Copy right hand side
      int m = matrix.RowLength;
      int n = matrix.ColumnLength;
      int nx = B.Columns;
      FloatMatrix ret = new FloatMatrix(n,nx);

#if MANAGED
      FloatMatrix X = new FloatMatrix(B);
      // Compute Y = transpose(Q)*B
      float[] column = new float[q_.RowLength];
      for (int j = 0; j < nx; j++) 
      {
        for (int k = 0; k < m; k++) 
        {
          column[k] = X.data[k][j];
        }
        for (int i = 0; i < m; i++) 
        {
          float s = 0;
          for (int k = 0; k < m; k++) 
          {
            s += q_.data[k][i] * column[k];
          }
          X.data[i][j] = s;
        } 
      }

      // Solve R*X = Y;
      for (int k = n-1; k >= 0; k--) 
      {
        for (int j = 0; j < nx; j++) 
        {
          X.data[k][j] /= r_.data[k][k];
        }
        for (int i = 0; i < k; i++) 
        {
          for (int j = 0; j < nx; j++) 
          {
            X.data[i][j] -= X.data[k][j]*r_.data[i][k];
          }
        }
      }
      for( int i = 0; i < n; i++ )
      {
        for( int j = 0; j < nx; j++ )
        {
          ret.data[i][j] = X.data[i][j];
        }
      }

#else
      float[] c = FloatMatrix.ToLinearArray(B);
      Lapack.Ormqr.Compute(Lapack.Side.Left, Lapack.Transpose.Trans, m, nx, n, qr, m, tau, c, m);
      Blas.Trsm.Compute(Blas.Order.ColumnMajor, Blas.Side.Left, Blas.UpLo.Upper, Blas.Transpose.NoTrans, Blas.Diag.NonUnit,
        n, nx, 1, qr, m, c, m);
      for ( int i = 0; i < n; i++ ) {
        for ( int j = 0; j < nx; j++) {
          ret.data[j*n+i] = c[j*m+(jpvt[i]-1)];
        }
      }

#endif
      return ret;
    }

 ///<summary>Constructor for QR decomposition class. The constructor performs the factorization and the upper and
 ///lower matrices are accessible by the <c>Q</c> and <c>R</c> properties.</summary>
 ///<param name="matrix">The matrix to factor.</param>
 ///<exception cref="ArgumentNullException">matrix is null.</exception>
 public FloatQRDecomp(IROFloatMatrix matrix) 
 {
   if (matrix == null)
     throw new System.ArgumentNullException("matrix cannot be null.");
   this.matrix = new FloatMatrix(matrix);
 }