C# LinearAlgebra.FloatVector类代码示例

    /// <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()
 {
   FloatVector test = new FloatVector(2);
   
   Assert.AreEqual(test.Length, 2);
   Assert.AreEqual(test[0],0);
   Assert.AreEqual(test[1],0);
 }

 public void CtorInitialValues()
 {
   FloatVector test = new FloatVector(2,1);
   
   Assert.AreEqual(test.Length, 2);
   Assert.AreEqual(test[0],1);
   Assert.AreEqual(test[1],1);
 }

 public void CtorArray()
 {
   float[] testvector = new float[2]{0,1};
   
   FloatVector test = new FloatVector(testvector);
   Assert.AreEqual(test.Length,testvector.Length);
   Assert.AreEqual(test[0],testvector[0]);
   Assert.AreEqual(test[1],testvector[1]);
 }

 public void CurrentException2() 
 {
   FloatVector test = new FloatVector(new float[2]{1f,2f});
   IEnumerator enumerator = test.GetEnumerator();
   enumerator.MoveNext();
   enumerator.MoveNext();
   enumerator.MoveNext();
   object value=enumerator.Current;
 }

 public void Current()
 {
   FloatVector test = new FloatVector(new float[2]{1f,2f});
   IEnumerator enumerator = test.GetEnumerator();
   bool movenextresult;
   
   movenextresult=enumerator.MoveNext();
   Assert.IsTrue(movenextresult);
   Assert.AreEqual(enumerator.Current,test[0]);
   
   movenextresult=enumerator.MoveNext();
   Assert.IsTrue(movenextresult);
   Assert.AreEqual(enumerator.Current,test[1]);
   
   movenextresult=enumerator.MoveNext();
   Assert.IsFalse(movenextresult);
 }

 public void SolveVector()
 {
   FloatVector b = new FloatVector(3);
   b[0] = 2;
   b[1] = 13;
   b[2] = 25;
   FloatVector x = cd.Solve(b);
   Assert.AreEqual(x[0],-3,TOLERENCE);
   Assert.AreEqual(x[1],8,TOLERENCE);
   Assert.AreEqual(x[2],8.333,TOLERENCE);
 }

		///<summary>Solves a system on linear equations, AX=B, where A is the factored matrixed.</summary>
		///<param name="B">RHS side of the system.</param>
		///<returns>the solution vector, X.</returns>
		///<exception cref="ArgumentNullException">B is null.</exception>
		///<exception cref="NotPositiveDefiniteException">A is not positive definite.</exception>
		///<exception cref="ArgumentException">The number of rows of A and the length of B must be the same.</exception>
		public FloatVector Solve(IROFloatVector B)
		{
			if (B == null)
			{
				throw new System.ArgumentNullException("B cannot be null.");
			}
			Compute();
			if (!ispd)
			{
				throw new NotPositiveDefiniteException();
			}
			else
			{
				if (B.Length != order)
				{
					throw new System.ArgumentException("The length of B must be the same as the order of the matrix.");
				}
#if MANAGED
				// Copy right hand side.
				FloatVector X = new FloatVector(B);
				// Solve L*Y = B;
				for (int i = 0; i < order; i++)
				{
					float sum = B[i];
					for (int k = i - 1; k >= 0; k--)
					{
						sum -= l.data[i][k] * X.data[k];
					}
					X.data[i] = sum / l.data[i][i];
				}
				// Solve L'*X = Y;
				for (int i = order - 1; i >= 0; i--)
				{
					float sum = X.data[i];
					for (int k = i + 1; k < order; k++)
					{
						sum -= l.data[k][i] * X.data[k];
					}
					X.data[i] = sum / l.data[i][i];
				}

				return X;
#else
                float[] rhs = FloatMatrix.ToLinearArray(B);
                Lapack.Potrs.Compute(Lapack.UpLo.Lower,order,1,l.data,order,rhs,B.Length);
                FloatVector ret = new FloatVector(order,B.Length);
                ret.data = rhs;
                return ret;
#endif
			}
		}

 ///<summary>Add a <c>FloatVector</c> to another <c>FloatVector</c></summary>
 ///<param name="lhs"><c>FloatVector</c> to add to.</param>
 ///<param name="rhs"><c>FloatVector</c> to add.</param>
 ///<returns><c>FloatVector</c> with results.</returns>
 public static FloatVector operator +(FloatVector lhs, FloatVector rhs)
 {
   FloatVector ret = new FloatVector(lhs);
   Blas.Axpy.Compute(ret.Length, 1, rhs.data, 1, ret.data, 1);
   return ret;
 }

    /// <summary>
    /// Solve the Yule-Walker equations for a symmetric square Toeplitz system
    /// </summary>
    /// <param name="R">The left-most column of the Toeplitz matrix.</param>
    /// <returns>The solution vector.</returns>
    /// <exception cref="ArgumentNullException">
    /// <B>R</B> is a null reference.
    /// </exception>
    /// <exception cref="ArgumentOutOfRangeException">
    /// The length of <B>R</B> must be greater than one.
    /// </exception>
    /// <exception cref="SingularMatrixException">
    /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
    /// </exception>
    /// <remarks>
    /// This member is used to solve the Yule-Walker system <B>AX</B> = -<B>a</B>,
    /// where <B>A</B> is a symmetric square Toeplitz matrix, constructed
    /// from the elements <B>R[0]</B>, ..., <B>R[N-2]</B> and
    /// the vector <B>a</B> is constructed from the elements
    /// <B>R[1]</B>, ..., <B>R[N-1]</B>.
    /// <para>
    /// Durbin's algorithm is used to solve the linear system. It requires
    /// approximately the <b>N</b> squared FLOPS to calculate the
    /// solution (<b>N</b> is the matrix order).
    /// </para>
    /// </remarks>
    public static FloatVector YuleWalker(IROFloatVector R)
    {

      FloatVector a;

      // check parameters
      if (R == null)
      {
        throw new System.ArgumentNullException("R");
      }
      else if (R.Length < 2)
      {
        throw new System.ArgumentOutOfRangeException("R", "The length of R must be greater than 1.");
      }
      else
      {

        int N = R.Length - 1;
        a = new FloatVector(N);                    // prediction coefficients
        FloatVector Z = new FloatVector(N);   // temporary storage vector
        float e;                            // predictor error
        float inner;                               // inner product
        float g;                                   // reflection coefficient
        int i, j, l;

        // setup first order solution
        e = R[0];
        if (e == 0.0f)
        {
          throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
        }
        g = -R[1] / R[0];
        a[0] = g;

        // calculate solution for successive orders
        for (i = 1; i < N; i++)
        {

          e *= (1.0f - g * g);
          if (e == 0.0f)
          {
            throw new SingularMatrixException("The Toeplitz matrix or one of the the leading sub-matrices is singular.");
          }

          // calculate inner product
          inner = R[i + 1];
          for (j = 0, l = i; j < i; j++, l--)
          {
            inner += a[j] * R[l];
          }

          // update prediction coefficients
          g = -(inner / e);
          for (j = 0, l = i - 1; j < i; j++, l--)
          {
            Z[j] = a[j] + g * a[l];
          }

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

          a[i] = g;

        }
      }

      return a;
    }

    /// <overloads>
    /// Solve a symmetric square Toeplitz system.
    /// </overloads>
    /// <summary>
    /// Solve a symmetric square Toeplitz system with a right-side vector.
    /// </summary>
    /// <param name="T">The left-most column of the Toeplitz matrix.</param>
    /// <param name="Y">The right-side vector of the system.</param>
    /// <returns>The solution vector.</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 length of <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
    /// vector and <B>Y</B> is a known vector.
    /// <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 FloatVector Solve(IROFloatVector T, IROFloatVector Y)
    {

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

        // allocate memory
        int N = T.Length;
        X = new FloatVector(N);                    // solution vector
        float e;                                   // prediction error

        // 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.");
        }
        X[0] = Y[0] / T[0];

        if (N > 1)
        {
          FloatVector a = new FloatVector(N - 1);   // prediction coefficients
          FloatVector Z = new FloatVector(N - 1);   // temporary storage vector
          float g;                                   // reflection coefficient
          float inner;                               // inner product
          float k;
          int i, j, l;

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

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

            a[i - 1] = g;

            e *= (1.0f - g * g);
            if (e == 0.0f)
            {
//.........这里部分代码省略.........

		///<summary>Implicit cast conversion to <c>DoubleVector</c> from <c>FloatVector</c></summary>
		static public DoubleVector ToDoubleVector(FloatVector src)
		{
			if (src == null)
			{
				return null;
			}
			DoubleVector ret = new DoubleVector(src.Length);
			Array.Copy(src.data, ret.data, src.Length);
			return ret;
		}

 ///<summary>Multiply a <c>FloatVector</c> x with a <c>float</c> y as x*y</summary>
 ///<param name="lhs"><c>FloatVector</c> as left hand operand.</param>
 ///<param name="rhs"><c>float</c> as right hand operand.</param>
 ///<returns><c>FloatVector</c> with results.</returns>
 public static FloatVector Multiply(FloatVector lhs, float rhs)
 {
   return lhs * rhs;
 }

 ///<summary>Add a <c>FloatVector</c> to another <c>FloatVector</c></summary>
 ///<param name="lhs"><c>FloatVector</c> to add to.</param>
 ///<param name="rhs"><c>FloatVector</c> to add.</param>
 ///<returns><c>FloatVector</c> with results.</returns>
 public static FloatVector Add(FloatVector lhs, FloatVector rhs)
 {
   return lhs + rhs;
 }

 ///<summary>Subtract a <c>FloatVector</c> from this<c>FloatVector</c></summary>
 ///<param name="vector"><c>FloatVector</c> to add.</param>
 ///<exception cref="ArgumentNullException">Exception thrown if null given as argument.</exception>
 public void Subtract(FloatVector vector)
 {
   if (vector == null)
   {
     throw new System.ArgumentNullException("FloatVector cannot be null.");
   }
   Blas.Axpy.Compute(this.Length, -1, vector.data, 1, this.data, 1);
 }