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)