matMulRectBand¶

Computes the transpose of a matrix, a matrix-vector product, or a matrix-matrix product, all matrices stored in band form.

Synopsis¶

matMulRectBand (string)

Required Arguments¶

char string (Input): String indicating matrix multiplication to be performed.

Return Value¶

The result of the multiplication is returned.

Optional Arguments¶

aMatrix, int nrowa, int ncola, int nlca, int nuca, float a (Input): The sparse matrix

\[A \in \Re^{\mathrm{nrowa} \times \mathrm{ncola}}\]

bMatrix, int nrowb, int ncolb, int nlcb, int nucb, float b (Input): The sparse matrix

\[B \in \Re^{\mathrm{nrowb} \times \mathrm{ncolb}}\]

xVector, int nx, float x, (Input): The vector x of length nx.
returnMatrixCodiagonals, nlc_result, nuc_result, (Output): If the function matMulRectBand returns data for a band matrix, use this option to retrieve the number of lower and upper codiagonals of the return matrix.

Description¶

The function matMulRectBand computes a matrix-matrix product or a matrix-vector product, where the matrices are specified in band format. The operation performed is specified by string. For example, if “A*x” is given, Ax is computed. In string, the matrices A and B and the vector x can be used. Any of these names can be used with trans, indicating transpose. The vector x is treated as a dense n × 1 matrix. If string contains only one item, such as “x” or “trans(A)”, then a copy of the array, or its transpose is returned.

The matrices and/or vector referred to in string must be given as optional arguments. Therefore, if string is “A*x”, then aMatrix and xVector must be given.

Examples¶

Example 1¶

Consider the matrix

\[\begin{split}A = \begin{bmatrix} 2 & -1 & 0 & 0 \\ -3 & 1 & -2 & 0 \\ 0 & 0 & -1 & 2 \\ 0 & 0 & 2 & 1 \\ \end{bmatrix}\end{split}\]

After storing A in band format, multiply A by \(x=(1,2,3,4)^T\) and print the result.

from pyimsl.math.matMulRectBand import matMulRectBand
from pyimsl.math.writeMatrix import writeMatrix

a = [[0., -1., -2., 2.], [2., 1., -1., 1.], [-3., 0., 2., 0.]]
n = 4
nlca = 1
nuca = 1
x = [1., 2., 3., 4.]

# Set b = A*x
b = matMulRectBand("A*x",
                   aMatrix={'nrowa': n, 'ncola': n,
                            'nlca': nlca, 'nuca': nuca, 'a': a},
                   xVector=x)

writeMatrix("Product, Ax", b)

Output¶

 
                    Product, Ax
          1            2            3            4
          0           -7            5           10

Example 2¶

This example uses the power method to determine the dominant eigenvector of E(100, 10). The same computation is performed by using eigSym, described in the chapter Eigensystem Analysis. The iteration stops when the component-wise absolute difference between the dominant eigenvector found by eigSym and the eigenvector at the current iteration is less than the square root of machine unit roundoff.

from __future__ import print_function
from numpy import *
from pyimsl.math.eigSym import eigSym
from pyimsl.math.generateTestBand import generateTestBand
from pyimsl.math.matMulRectBand import matMulRectBand
from pyimsl.math.vectorNorm import vectorNorm
from pyimsl.math.machine import machine
from math import sqrt

n = 100
c = 10
tolerance = sqrt(machine(4))
error = 1.0
evec = zeros((n), dtype='double')
z = zeros((n), dtype='double')
q = zeros((n), dtype='double')
dense_a = zeros((n, n), dtype='double')
a = generateTestBand(n, c)

# Convert to dense format, starting with upper triangle
start = c
for i in range(0, c):
    k = 0
    for j in range(start, n):
        dense_a[k, j] = a[i, j]
        k = k + 1
    start = start - 1

# Convert diagonal
for j in range(0, n):
    dense_a[j, j] = a[c, j]

# Convert lower triangle
stop = n - 1
for i in range(c + 1, 2 * c + 1):
    k = i - c
    for j in range(0, stop):
        dense_a[k, j] = a[i, j]
        k = k + 1
    stop = stop - 1

# Determine dominant eigenvector by a dense method
dense_evec = []
dense_eval = eigSym(dense_a, vectors=dense_evec)
for i in range(0, n):
    evec[i] = dense_evec[i][0]

# Normalize
norm = vectorNorm(evec)
for i in range(0, n):
    evec[i] = evec[i] / norm
for i in range(0, n):
    q[i] = 1.0 / sqrt(float(n))

# Do power method
index = []
while (error > tolerance):
    z = matMulRectBand("A*x",
                       aMatrix={'nrowa': n, 'ncola': n,
                                'nlca': c, 'nuca': c, 'a': a},
                       xVector=q)

    # Normalize
    z2 = z.reshape(n)
    norm = vectorNorm(z2)
    for i in range(0, n):
        q[i] = z2[i] / norm

    # Compute maximum absolute error between
    # any two elements
    error = vectorNorm(q, secondVector=evec, infNorm=index)
print("Maximum absolute error = ", error)

Output¶

Maximum absolute error =  1.4413304612326261e-08