Triangular Matrices in IronPython QuickStart Sample

Illustrates how to work efficiently with upper or lower triangular or trapezoidal matrices in IronPython.

This sample is also available in: C#, Visual Basic, F#.

Overview

This QuickStart sample demonstrates how to work with triangular matrices in Numerics.NET. Triangular matrices are special matrices where all elements either above the diagonal (lower triangular) or below the diagonal (upper triangular) are zero.

The sample covers:

Creating upper and lower triangular matrices through various constructors
Working with unit diagonal triangular matrices (matrices with all 1’s on the diagonal)
Extracting triangular parts from dense matrices
Accessing matrix properties like IsLowerTriangular and IsUpperTriangular
Handling matrix elements and the restrictions on modifying zero elements
Working with rows and columns of triangular matrices

The code illustrates proper error handling when attempting to modify read-only elements and shows how to efficiently work with these specialized matrix types. This is particularly useful in numerical algorithms where triangular matrices arise naturally, such as in LU decomposition or backward substitution.

The code

from System import Array

import numerics

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

# Illustrates the use of the TriangularMatrix class in the 
# Extreme.Mathematics.LinearAlgebra namespace of Extreme Numerics.NET.

# Triangular matrices are matrices whose elements
# above or below the diagonal are all zero. The
# former is called lower triangular, the latter
# lower triangular. In addition, triangular matrices
# can have all 1's on the diagonal.

#
# Constructing triangular matrices
#

# Constructing triangular matrices is similar to
# constructing general matrices. See the
# BasicMatrices QuickStart samples for a more
# complete discussion.
#
# All constructors take a MatrixTriangle
# value as their first parameter. This indicates
# whether an upper or lower triangular matrix 
# should be created. The following creates a
# 5x5 lower triangular matrix:
t1 = Matrix.CreateLowerTriangular(5, 5)
# You can also specify whether the diagonal 
# consists of all 1's using a unitDiagonal parameter:
t2 = Matrix.CreateLowerTriangular(5, 5, MatrixDiagonal.UnitDiagonal)
# Triangular matrices access and modify only the
# elements that are non-zero. If the diagonal
# mode is UnitDiagonal, the diagonal elements
# are not used, since they are all equal to 1.
components = Array[float]([ \
    11, 12, 13, 14, 15, \
    21, 22, 23, 24, 25, \
    31, 32, 33, 34, 35, \
    41, 42, 43, 44, 45, \
    51, 52, 53, 54, 55 ])
# The following creates a matrix using the
# upper triangular part of the above.
t3 = Matrix.CreateUpperTriangular(5, 5, components, MatrixElementOrder.RowMajor)
print "t3 =", t3
# Same as above, but unit diagonal:
t4 = Matrix.CreateUpperTriangular(5, 5, components, MatrixDiagonal.UnitDiagonal, MatrixElementOrder.RowMajor, True)
print "t4 =", t4

#
# Extracting triangular matrices
#

# You may want to use part of a dense matrix
# as a triangular matrix. The static 
# ExtractUpperTriangle and ExtractLowerTriangle
# methods perform this task.
m = Matrix.Create(5, 5, components, MatrixElementOrder.ColumnMajor)
print "m =", m
# Both methods are overloaded. The simplest
# returns a triangular matrix of the same dimension:
t5 = TriangularMatrix.ExtractLowerTriangle(m)
print "t5 =", t5
# You can also specify if the matrix is unit diagonal:
t6 = TriangularMatrix.ExtractUpperTriangle(m, MatrixDiagonal.UnitDiagonal)
print "t6 =", t6
# Or the dimensions of the matrix if they don't 
# match the original:
t7 = TriangularMatrix.ExtractUpperTriangle(m, 3, 3, MatrixDiagonal.UnitDiagonal)
print "t7 =", t7
print 

#
# TriangularMatrix properties
#

# The IsLowerTriangular and IsUpperTriangular return
# a boolean value:
print "t4 is lower triangular? -", t4.IsLowerTriangular
print "t4 is upper triangular? -", t4.IsUpperTriangular
# The IsUnitDiagonal property indicates whether the
# matrix has all 1's on its diagonal:
print "t3 is unit diagonal? -", t3.IsUnitDiagonal
print "t4 is unit diagonal? -", t4.IsUnitDiagonal
print 
# You can get and set matrix elements:
t3[1, 3] = 55
print "t3[1, 3] =", t3[1, 3]
# But trying to set an element that is zero or
# is on the diagonal for a unit diagonal matrix
# causes an exception to be thrown:
try:
	t3[3, 1] = 100
except ComponentReadOnlyException as e:
	print "Error accessing element:", e.Message

#
# Rows and columns
#

# The GetRow and GetColumn methods are
# available.
row = t3.GetRow(1)
row = t3[1,:]
print "row 2 of t3 =", row
column = t4.GetColumn(1, 0, 2)
column = t4[0:3, 1]
print "2nd column of t4 from row 1 to 3 =", column