matMulRectCoordinate¶

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

Synopsis¶

matMulRectCoordinate (string)

Required Arguments¶

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

Return Value¶

The returned value is the result of the multiplication. If the result is a vector, the return type is float. If the result of the multiplication is a sparse matrix, the return type is sparse_elem.

Optional Arguments¶

aMatrix, int nrowa, int ncola, int nza, sparse_elem a (Input)

The sparse matrix

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

with nza nonzero elements.

bMatrix, int nrowb, int ncolb, int nzb, sparse_elem b (Input)

The sparse matrix

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

with nzb nonzero elements.

xVector, float x, (Input)

The vector x of length nx.

returnMatrixSize (Output)

If the function matMulRectCoordinate returns a vector of type sparse_elem, use this option to retrieve the length of the return vector, i.e. the number of nonzero elements in the sparse matrix generated by the requested computations.

sparseOutput (Output)

If the result of the computation in a vector, return the answer in the user supplied sparse sparseOutput. It’s size depends on the computation.

Description¶

The function matMulRectCoordinate computes a matrix-matrix product or a matrix-vector product, where the matrices are specified in coordinate representation. 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. Some multiplications, such as “A*trans(A)” or “trans(x)*B”, will produce a sparse matrix in coordinate format as a result. Other products such as “B*x” will produce a floating type, containing the resulting vector.

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¶

In this example, a sparse matrix in coordinate form is multipled by a vector.

from numpy import *
from pyimsl.math.matMulRectCoordinate import matMulRectCoordinate
from pyimsl.math.writeMatrix import writeMatrix

a = [[0, 0, 10.0],
     [1, 1, 10.0],
     [1, 2, -3.0],
     [1, 3, -1.0],
     [2, 2, 15.0],
     [3, 0, -2.0],
     [3, 3, 10.0],
     [3, 4, -1.0],
     [4, 0, -1.0],
     [4, 3, -5.0],
     [4, 4, 1.0],
     [4, 5, -3.0],
     [5, 0, -1.0],
     [5, 1, -2.0],
     [5, 5, 6.0]]
n = 6
b = array([10.0, 7.0, 45.0, 33.0, -34.0, 31.0])

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

writeMatrix("Product Ab", x)

Output¶

 
                                 Product Ab
          1            2            3            4            5            6
        100          -98          675          344         -302          162

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. 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.generateTestCoordinate import generateTestCoordinate
from pyimsl.math.matMulRectCoordinate import matMulRectCoordinate
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 = generateTestCoordinate(n, c)

# Convert to dense format
for i in range(0, len(a)):
    dense_a[a[i][0], a[i][1]] = a[i][2]

# 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 = matMulRectCoordinate("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.441330455681511e-08