13,861,429 members

#### Stats

83.8K views
92 bookmarked
Posted 19 Sep 2012
Licenced CPOL

# Yet Another Math Parser (YAMP)

, 30 Sep 2012
Constructing a fast math parser using Reflection to do numerics like Matlab.
 ```using System; using YAMP; namespace YAMP.Numerics { /// /// Cholesky Decomposition. /// For a symmetric, positive definite matrix A, the Cholesky decomposition /// is an lower triangular matrix L so that A = L*L'. /// If the matrix is not symmetric or positive definite, the constructor /// returns a partial decomposition and sets an internal flag that may /// be queried by the isSPD() method. /// public class CholeskyDecomposition : DirectSolver { #region Members /// /// Array for internal storage of decomposition. /// double[][] L; /// /// Row and column dimension (square matrix). /// int n; /// /// Symmetric and positive definite flag. /// bool isspd; #endregion // Class variables #region Constructor /// /// Cholesky algorithm for symmetric and positive definite matrix. /// /// Square, symmetric matrix. /// Structure to access L and isspd flag. public CholeskyDecomposition(MatrixValue Arg) { // Initialize. double[][] A = Arg.GetRealArray(); n = Arg.DimensionY; L = new double[n][]; for (int i = 0; i < n; i++) L[i] = new double[n]; isspd = (Arg.DimensionX == n); // Main loop. for (int j = 0; j < n; j++) { double[] Lrowj = L[j]; double d = 0.0; for (int k = 0; k < j; k++) { double[] Lrowk = L[k]; double s = 0.0; for (int i = 0; i < k; i++) s += Lrowk[i] * Lrowj[i]; Lrowj[k] = s = (A[j][k] - s) / L[k][k]; d = d + s * s; isspd = isspd & (A[k][j] == A[j][k]); } d = A[j][j] - d; isspd = isspd & (d > 0.0); L[j][j] = Math.Sqrt(Math.Max(d, 0.0)); for (int k = j + 1; k < n; k++) L[j][k] = 0.0; } } #endregion // Constructor #region Properties /// /// Is the matrix symmetric and positive definite? /// /// true if A is symmetric and positive definite. public bool SPD { get { return isspd; } } #endregion // Public Properties #region Public Methods /// /// Return triangular factor. /// /// L public virtual MatrixValue GetL() { return new MatrixValue(L, n, n); } /// Solve A*X = B /// A Matrix with as many rows as A and any number of columns. /// /// X so that L*L'*X = B /// /// Matrix row dimensions must agree. /// /// Matrix is not symmetric positive definite. /// public override MatrixValue Solve(MatrixValue B) { if (B.DimensionY != n) throw new DimensionException(n, B.DimensionY); if (!isspd) throw new MatrixFormatException("symmetric positive definite"); // Copy right hand side. double[][] X = B.GetRealArray(); int nx = B.DimensionX; // Solve L*Y = B; for (int k = 0; k < n; k++) { for (int i = k + 1; i < n; i++) { for (int j = 0; j < nx; j++) X[i][j] -= X[k][j] * L[i][k]; } for (int j = 0; j < nx; j++) X[k][j] /= L[k][k]; } // Solve L'*X = Y; for (int k = n - 1; k >= 0; k--) { for (int j = 0; j < nx; j++) X[k][j] /= L[k][k]; for (int i = 0; i < k; i++) { for (int j = 0; j < nx; j++) X[i][j] -= X[k][j] * L[k][i]; } } return new MatrixValue(X, n, nx); } #endregion // Public Methods } }```

If a file you wish to view isn't highlighted, and is a text file (not binary), please let us know and we'll add colourisation support for it.

## Share

 Architect Germany
Florian lives in Munich, Germany. He started his programming career with Perl. After programming C/C++ for some years he discovered his favorite programming language C#. He did work at Siemens as a programmer until he decided to study Physics.

During his studies he worked as an IT consultant for various companies. After graduating with a PhD in theoretical particle Physics he is working as a senior technical consultant in the field of home automation and IoT.

Florian has been giving lectures in C#, HTML5 with CSS3 and JavaScript, software design, and other topics. He is regularly giving talks at user groups, conferences, and companies. He is actively contributing to open-source projects. Florian is the maintainer of AngleSharp, a completely managed browser engine.