Structured Linear Equations in IronPython QuickStart Sample
Illustrates how to solve systems of simultaneous linear equations that have special structure in IronPython.
View this sample in: C# Visual Basic F#
```Python import numerics from System import Array from Extreme.Mathematics import * # The structured matrix classes reside in the # Extreme.Mathematics.LinearAlgebra namespace. from Extreme.Mathematics.LinearAlgebra import * # Illustrates solving symmetrical and triangular systems # of simultaneous linear equations using classes # in the Extreme.Mathematics.LinearAlgebra namespace of the Extreme # Optimization Mathematics Library for .NET. # To learn more about solving general systems of # simultaneous linear equations, see the # LinearEquations QuickStart Sample. # # The methods and classes available for solving # structured systems of equations are similar # to those for general equations. # # Triangular systems and matrices # print "Triangular matrices:" # For the basics of working with triangular # matrices, see the TriangularMatrices QuickStart # Sample. # # Let's start with a triangular matrix. Remember # that elements are stored in column-major order # by default. t = Matrix.CreateUpperTriangular(4, 4, Array[float]([ \ 1, 0, 0, 0, \ 1, 2, 0, 0, \ 1, 4, 1, 0, \ 1, 3, 1, 2 ]), \ MatrixElementOrder.ColumnMajor) b1 = Vector([ 1, 3, 6, 3 ]) b2 = Matrix([[1,2], [3,3], [6,5], [3,8]]) print "t = {0:.4f}".format(t) # # The Solve method # # The following solves m x = b1. The second # parameter specifies whether to overwrite the # right-hand side with the result. x1 = t.Solve(b1, False) print "x1 = {0:.4f}".format(x1) # If the overwrite parameter is omitted, the # right-hand-side is overwritten with the solution: t.Solve(b1) print "b1 = {0:.4f}".format(b1) # You can solve for multiple right hand side # vectors by passing them in a DenseMatrix: x2 = t.Solve(b2, False) print "x2 = {0:.4f}".format(x2) # # Related Methods # # You can verify whether a matrix is singular # using the IsSingular method: print "IsSingular(t) =", t.IsSingular() # The inverse matrix is returned by the Inverse # method: print "Inverse(t) = {0:.4f}".format(t.GetInverse()) # The determinant is also available: print "Det(t) = {0:.4f}".format(t.GetDeterminant()) # The condition number is an estimate for the # loss of precision in solving the equations print "Cond(t) = {0:.4f}".format(t.EstimateConditionNumber()) print # # Symmetric systems and matrices # print "Symmetric matrices:" # For the basics of working with symmetric # matrices, see the SymmetricMatrices QuickStart # Sample. # # Let's start with a symmetric matrix. Remember # that elements are stored in column-major order # by default. s = Matrix.CreateSymmetric(4, Array[float]([ \ 1, 0, 0, 0, \ 1, 2, 0, 0, \ 1, 1, 2, 0, \ 1, 0, 1, 4 ]), MatrixTriangle.Upper) b1 = Vector.Create(1, 3, 6, 3) print "s = {0:.4f}".format(s) # # The Solve method # # The following solves m x = b1. The second # parameter specifies whether to overwrite the # right-hand side with the result. x1 = s.Solve(b1, False) print "x1 = {0:.4f}".format(x1) # If the overwrite parameter is omitted, the # right-hand-side is overwritten with the solution: s.Solve(b1) print "b1 = {0:.4f}".format(b1) # You can solve for multiple right hand side # vectors by passing them in a DenseMatrix: x2 = s.Solve(b2, False) print "x2 = {0:.4f}".format(x2) # # Related Methods # # You can verify whether a matrix is singular # using the IsSingular method: print "IsSingular(s) =", s.IsSingular() # The inverse matrix is returned by the Inverse # method: print "Inverse(s) = {0:.4f}".format(s.GetInverse()) # The determinant is also available: print "Det(s) = {0:.4f}".format(s.GetDeterminant()) # The condition number is an estimate for the # loss of precision in solving the equations print "Cond(s) = {0:.4f}".format(s.EstimateConditionNumber()) print # # The CholeskyDecomposition class # # If the symmetric matrix is positive definite, # you can use the CholeskyDecomposition class # to optimize performance if multiple operations # need to be performed. This class does the # bulk of the calculations only once. This # decomposes the matrix into G x transpose(G) # where G is a lower triangular matrix. # # If the matrix is indefinite, you need to use # the LUDecomposition class instead. See the # LinearEquations QuickStart Sample for details. print "Using Cholesky Decomposition:" # The constructor takes an optional second argument # indicating whether to overwrite the original # matrix with its decomposition: cf = s.GetCholeskyDecomposition(False) # The Factorize method performs the actual # factorization. It is called automatically # if needed. cf.Decompose() # All methods mentioned earlier are still available: x2 = cf.Solve(b2, False) print "x2 = {0:.4f}".format(x2) print "IsSingular(m) =", cf.IsSingular() print "Inverse(m) = {0:.4f}".format(cf.GetInverse()) print "Det(m) = {0:.4f}".format(cf.GetDeterminant()) print "Cond(m) = {0:.4f}".format(cf.EstimateConditionNumber()) # In addition, you have access to the # triangular matrix, G, of the composition. print " G = {0:.4f}".format(cf.LowerTriangularFactor) # Note that if the matrix is indefinite, # the factorization will fail and throw a # MatrixNotPositiveDefiniteException. s[0, 0] = -99 cf = s.GetCholeskyDecomposition() try: cf.Decompose() except MatrixNotPositiveDefiniteException as e: print e.Message # To avoid this, you can use the TryDecompose method, # which returns True if the decomposition succeeds, # and false otherwise: if cf.TryDecompose(): print "Decomposition succeeded!" else: print "Decomposition failed!" ```