Band Matrices QuickStart Sample (IronPython)

import numerics

# The BandMatrix class resides in the Extreme.Mathematics.LinearAlgebra
# namespace.
from Extreme.Mathematics import *
from Extreme.Mathematics.LinearAlgebra import *

#/ Illustrates the use of the BandMatrix class in the 
#/ Extreme.Mathematics.LinearAlgebra namespace of the Extreme Optimization
#/ Numerical Libraries for .NET.

# Band matrices are matrices whose elements
# are nonzero only in a diagonal band around
# the main diagonal.
#
# General band matrices, upper and lower band
# matrices, and symmetric band matrices are all
# represented by a single class: BandMatrix.

#
# Constructing band matrices
#

# Constructing band matrices is similar to
# constructing general matrices. It is done by
# calling a factory method on the Matrix class.
# See theBasicMatrices QuickStart samples 
# for a more complete discussion.

# The following creates a 7x5 band matrix with
# upper bandwidth 1 and lower bandwidth 2:
b1 = Matrix.CreateBanded(7, 5, 2, 1)

# Once the upper and lower bandwidth are set, # it cannot be changed. Elements that are outside
# the band cannot be set.

# A second factory method lets you create upper
# or lower band matrices. The following constructs
# an 11x11 upper band matrix with unit diagonal
# and three non-zero upper diagonals.
b2 = Matrix.CreateUpperBanded(11, 11, 3, MatrixDiagonal.UnitDiagonal)

# To create a symmetric band matrix, you only need
# the size and the bandwith. The following creates
# a 6x6 symmetric tri-diagonal matrix:
b3 = Matrix.CreateSymmetricBanded(7, 1)

# We can assign values to the components by using
# the GetDiagonal method.
b3.GetDiagonal(0).SetValue(2)
b3.GetDiagonal(1).SetValue(-1)

# Extracting band matrices

# Another way to construct a band matrix is by
# extracting them from an existing matrix.
m = Matrix([[1,3,5], [2,4,4], [3,3,5], [2,4,7]])
# To get the lower band part of m with bandwidth 2:
b4 = BandMatrix.Extract(m, 2, 0)

#
# BandMatrix properties
#

# A number of properties are available to determine 
# whether a BandMatrix has a special structure:
print "b2 is upper?", b2.IsUpperTriangular
print "b2 is lower?", b2.IsUpperTriangular
print "b2 is unit diagonal?", b2.IsUnitDiagonal
print "b2 is symmetrical?", b2.IsSymmetrical

#
# BandMatrix methods
#

# You can get and set matrix elements:
b3[2, 3] = 55
print "b3[2, 3] =", b3[2, 3]
# And the change will automatically be reflected
# in the symmetric element:
print "b3[3, 2] =", b3[3, 2]

#
# Row and column views
#

# The GetRow and GetColumn methods are
# available.
row = b2.GetRow(1)
row = b2[1,:]
print "row 1 of b2 =", row
column = b2.GetColumn(2, 3, 4)
column = b2[3:5,2]
print "column 3 of b2 from row 4 to row 5 =", column

#
# Band matrix decompositions
#

# Specialized classes exist to represent the
# LU decomposition of a general band matrix
# and the Cholesky decomposition of a 
# symmetric band matrix.

# Because of pivoting, the upper band matrix of
# the LU decomposition has larger bandwidth.
# You need to allocate extra space to be able to
# overwrite a matrix with its LU decomposition.

# The following creates a 7x5 band matrix with
# upper bandwidth 1 and lower bandwidth 2.
b5 = Matrix.CreateBanded(7, 7, 2, 1, True)
b5.GetDiagonal(0).SetValue(2.0)
b5.GetDiagonal(-2).SetValue(-1.0)
b5.GetDiagonal(1).SetValue(-1.0)

# Other than that, the API is the same as
# other decomposition classes.
blu = b5.GetLUDecomposition(True)
solution = blu.Solve(Vector.CreateConstant(b5.ColumnCount, 1.0))
print "solution of b5*x = ones: {0:.4f}".format(solution)