Structured Linear Equations QuickStart Sample (IronPython)

Illustrates how to solve systems of simultaneous linear equations that have special structure in IronPython.

C# code Visual Basic code F# code Back to QuickStart Samples

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!"