Jacob C. Kesinger > Math-MatrixSparse-0.01 > Math::MatrixSparse
Module Version: 0.01

# NAME

Math::MatrixSparse - Perl extension for sparse matrices.

# SYNOPSIS

Math::MatrixSparse is a module implementing naive sparse matrices. Its syntax is designed to partially overlap with Math::MatrixReal where possible and reasonable.

Basic matrix operations are present, including addition, subtraction, scalar multiplication, matrix multiplication, transposition, and exponentiation.

Three methods of solving systems iteratively are available, including Jacobi, Gauss-Seidel, and Symmetric Over-Relaxation.

Real-valued matrices can be read from files in the Matrix Market and Harwell Boeing formats, and written in the Matrix Market format. In addition, they can be read from a given string.

Certain special types of matrices are understood, but not optimized for, such as upper and lower triangular, diagonal, symmetric, skew-symmetric, positive, negative, and pattern. Methods are available to determine the properties of a given matrix.

Individual rows, columns, and diagonals of matrices can be obtained, as can the upper and lower triangular, symmetric, skew symmetric, positive and negative portions.

# DESCRIPTION

• `use Math::MatrixSparse;`

Load the module and make its methods and operators available.

## CREATION AND INPUT-OUTPUT METHODS

• `Math::MatrixSparse->new(\$name)`

`new Math::MatrixSparse(\$name)`

Creates a new empty matrix named \$name, which may be undef.

• `Math::MatrixSparse->newfromstring(\$string,\$name)`

Creates a new matrix named \$name based on \$string. \$string is of the pattern /(\d+)\s+(\d+)\s+(.+)\n/, where (\$row, \$column, \$value) = (\$1,\$2,\$3).

Example:

``` \$matrixspec = <<MATRIX;
1 1 1
1 2 2
2 1 3
3 3 4
MATRIX
my \$MS = Math::MatrixSparse->newfromstring(\$matrixspec,"MS");```

\$MS has four elements, (1,1), (1,2), (2,1), and (3,3), with the values 1,2,3, and 4.

• `Math::MatrixSparse->newdiag(\$arrayref,\$name)`

Creates a new square matrix named \$name from the elements in \$arrayref. \$arrayref[0] will become \$matrix(1,1), etc.

• `Math::MatrixSparse->newdiagfromstring(\$string, \$name)`

Similar to `Math::MatrixSparse->newfromstring` except that \$string is of the pattern /(\d+)\s+(.+)\n/, and (\$row, \$column, \$value) = (\$1, \$1, \$2).

Example:

``` \$diagspec = <<MATRIX;
1   1
2   -1
20  1
300 -1
MATRIX
my \$MD = Math::MatrixSparse->newdiagfromstring(\$matrixspec,"MS");```

\$MD has four elements, (1,1), (2,2), (20,20), and (300,300), with the values 1, -1, 1, and -1, and is square.

• `Math::MatrixSparse->newrandom(\$maxrow, \$maxcol, \$maxentries, \$density,\$name)`

Creates a new matrix of the specified size (\$maxrow*\$maxcol) with at most \$maxentries. If \$density is between 0 and 1 (inclusive), the expected number of entries is \$maxentries*\$density. If \$maxentries is undef, it is set to \$maxrow*\$maxcol; if \$maxcol is missing, it is set to \$maxrow. If \$density is missing, or outside the valid range, it is set to one.

• `Math::MatrixSparse->newmatrixmarket(\$filename)`

Creates a new matrix based on data stored in the file \$filename, in Matrix Market format. See http://www.nist.gov/MatrixMarket/ for more information on this format.

Pattern matrix data has elements set to one. Matrices flagged as symmetric or skew symmetric have symmetrify() or skewsymmetrify() applied.

Exits with an error if \$filename cannot be read from.

• `Math::MatrixSparse->newharwellboeing(\$filename)`

Creates a new matrix based on data stored in the file \$filename, in Harwell-Boeing format. See http://math.nist.gov/MatrixMarket/collections/hb.html for more information on this format.

Pattern matrix data has elements set to one. Matrices flagged as symmetric or skew symmetric have symmetrify() or skewsymmetrify() applied.

Exits with an error if \$filename cannot be read from.

• `\$matrix->writematrixmarket(\$filename)`

The contents of \$matrix are written, in Matrix Market format, to \$filename. It is assumed that \$matrix->{rows} and \$matrix->{columns} are accurate--\$matrix->sizetofit() should be called if this is not the case. No optimizations are performed for symmetric, skew symmetric, or pattern data, and real values are assumed.

See http://www.nist.gov/MatrixMarket/ for more information on this format.

Exits with an error if \$filename cannot be written to.

• `Math::MatrixSparse->newidentity(\$n,\$m)`

Creates an identity matrix with \$n rows and \$m columns. If \$m is omitted it is assumed equal to \$n.

• `print \$matrix`

`\$matrix->print(\$name)`

`\$matrix->print()`

Prints the element in the matrix in the format \$name(\$row,\$colum) = \$value

Calling this as a method with an argument overrides the value in \$matrix->{name}.

## INTERFACE METHODS

• `(\$i,\$j) = &splitkey(\$key)`

`\$key = &makekey(\$i,\$j)`

These two routines convert a position (\$i,\$j) of the matrix into a key \$key. Programs that use Math::MatrixSparse should not have to use keys at all, but if they do, they should use these routines.

## GENERAL METHODS

• `\$matrix2 = \$matrix1->copy()`

`\$matrix2 = \$matrix1->copy(\$name)`

Returns an exact copy of \$matrix1, including dimensions, name, and special flags. If \$name is given, it will be the name of \$matrix2.

• `\$matrix->name(\$name)`

Gives a name to \$matrix. Useful mostly in printing. Certain methods attempt to create useful names, but you probably want to set your own.

• `\$matrix->assign(\$i,\$j,\$value)`

Puts \$value in row \$i, column \$j of \$matrix. Updates \$matrix->{row} and \$matrix->{column} if necessary. Removes or modifies special flags when necessary.

• `\$matrix->assignspecial(\$i,\$j,\$value)`

As \$matrix->assign(), except that it preserves symmetry and skew-symmetry. That is, if \$matrix is marked symmetric, assigning a value to (\$i,\$j) also assigns the value to (\$j,\$i). Also preserves patterns.

• `\$matrix->assignkey(\$key,\$value)`

As \$matrix->assign(), except \$i and \$j have already been combined into \$key.

• `\$matrix->element(\$i,\$j)`

Returns the element at row \$i column \$j of \$matrix.

• `\$matrix->elementkey(\$key)`

As \$matrix->element(), except \$i and \$j have already been combined into \$key.

• `\$matrix->elements()`

Returns an array consisting of all the elements of \$matrix, in key form, suitable for a loop.

See also SORTING METHODS below if the order of the elements is important.

Example: foreach my \$key (\$matrix->elements()) { my (\$i,\$j) = &splitkey(\$key); ... }

• `\$matrix->delete(\$i,\$j)`

Deletes the (\$i,\$j) element of \$matrix.

• `\$matrix->deletekey(\$key)`

As \$matrix->delete(), except \$i and \$j have already been combined into \$key.

• `\$matrix->cull()`

Returns a copy of \$matrix with any zero elements removed. Equivalent to \$matrix->threshold(0).

• `\$matrix->threshold(\$thresh)`

Returns a copy of \$matrix with any elements with absolute value less than \$thresh removed.

• `\$matrix->sizetofit()`

Returns a copy of \$matrix with dimensions sufficient only to cover its data.

• `\$matrix->clearspecials()`

Removes all knowledge of special properties of \$matrix. Use \$matrix->diagnose() to put them back.

## DECOMPOSITIONAL METHODS

• `\$matrix->row(\$i)`

`\$matrix->row(\$i,\$persist)`

Returns a matrix consisting only of the \$i th row of \$matrix. If \$persist is non-zero, the elements of \$matrix are sorted by row, to speed up future calls to row().

Note that many methods, most notably assign() and delete(), remove the data necessary to use the fast algorithm. If this data is present, a binary search is used instead of an exhaustive term-by-term search. There is an intial cost of O(n log n) operations (n==\$matrix->count()) to use this data, but each search costs O(log n) operations. This is in comparison to O(n) operations without the data.

So, in summary, if many separate rows are needed, and the matrix is unchanging, the first (at least) call to row() should have a non-zero value of \$persist.

• `\$matrix->column(\$j)`

`\$matrix->column(\$j,\$persist)`

Returns a matrix consisting only of the \$j th column of \$matrix. See also \$matrix->row().

• `\$matrix->diagpart(\$offset)`

If \$offset is zero or undefined, returns a matrix consisting of the diagonal elements of \$matrix, if any. If \$offset is non-zero, returns a parallel diagonal. For example, if \$offset=1, the first superdiagonal is returned. Posive \$offset s return diagonals above the main diagonal, negative offsets return diagonals below the main diagonal.

The returned matrix is the same size as \$matrix.

• `\$matrix->lowerpart()`

Returns a matrix consisting of the strictly lower triangular parts of \$matrix, if any. The returned matrix is the same size as \$matrix.

• `\$matrix->upperpart()`

Returns a matrix consisting of the strictly upper triangular parts of \$matrix, if any. The returned matrix is the same size as \$matrix.

• `\$matrix->nondiagpart(\$offset)` `\$matrix->nonlowerpart()` `\$matrix->nonupperpart()`

As before, except that the returned matrix is everything except the part specified.

Example:

``` A =  D1  U1  U2
L1  D2  U3
L2  L3  D3

\$A->lowerpart() == L1 L2 L3
\$A->diagpart() == D1 D2 D3
\$A->diagpart(1) == U1 U3
\$A->diagpart(-1) == L1 L3
\$A->upperpart() == U1 U2 U3
\$A->nonlowerpart() == D1 U1 U2 D2 U3 D3
\$A->nondiagpart(1) == D1 U2 L1 D2 L2 L3 D3```
• `\$matrix->symmetricpart()`

Returns the symmetric part of \$matrix, i.e. 0.5*(\$matrix+\$matrix->transpose()).

• `\$matrix->skewsymmetricpart()`

Returns the skewsymmetric part of \$matrix, i.e. 0.5*(\$matrix-\$matrix->transpose()).

• `\$matrix->positivepart()`

Returns the positive part of \$matrix, i.e. those elements of \$matrix larger than zero.

• `\$matrix->negativepart()`

Returns the positive part of \$matrix, i.e. those elements of \$matrix larger than zero.

• `\$matrix->mask(\$i1,\$i2,\$j1,\$j2)`

Returns a matrix whose only elements are inside (\$i1,\$j1) (\$i1,\$j2) (\$i2,\$j1) (\$i2,\$j2)

• `\$matrix->submatrix(\$i1,\$i2,\$j1,\$j2)`

Returns that portion of \$matrix between (\$i1,\$j1) (\$i1,\$j2) (\$i2,\$j1) (\$i2,\$j2)

Subscripts are changed, so \$matrix(\$i1,\$j1) is \$submatrix(1,1).

• `\$matrix1->insert(\$i,\$j,\$matrix2)`

The values of \$matrix2 are assigned to the values of \$matrix1, offset by (\$i,\$j). That is, \$matrix1(\$i+\$k,\$j+\$l)=\$matrix2(\$j,\$l).

• `\$A = \$matrix->shift(\$row,\$col)`

Creates a new matrix \$A, \$matrix(\$i,\$j) == \$A(\$i+\$row,\$j+\$col)

## INFORMATIONAL METHODS

• `\$matrix->dim()`

Returns (\$matrix->{rows}, \$matrix->{columns})

• `\$matrix->exactdim()`

Calculates explicitly the dimension of \$matrix. This differs from \$matrix->dim() in that \$matrix->{rows} and \$matrix->{columns} may not have been updated properly.

• `\$matrix->count()`

Returns the number of defined elements in \$matrix, 0 or otherwise.

• `\$matrix->width()`

Calculates the bandwidth of \$matrix, i.e. the maximum value of abs(\$i-\$j) for all elements at (\$i,\$j) in \$matrix.

• `\$matrix->sum()`

Finds the sum of all elements in \$matrix.

• `\$matrix->abssum()`

Finds the sum of the absolute values of all elements in \$matrix.

• `\$matrix->sumeach(\$coderef)`

Applies &\$coderef to each of the elements of \$matrix, and sums the result. See `\$matrix->each(\$coderef)` for more details.

• `\$matrix->columnnorm()`

`\$matrix->norm_one()`

Finds the maximum column sum of \$matrix. That is, for each column of \$matrix, find the sum of absolute values of its elements, then return the largest such sum. `norm_one` is provided for compatibility with Math::MatrixReal.

• `\$matrix->rownorm()`

`\$matrix->norm_max()`

Finds the maximum row sum of \$matrix. That is, for each row of \$matrix, find the sum of absolute values of its elements, then return the largest such sum. `norm_max` is provided for compatibility with Math::MatrixReal.

• `\$matrix->norm_frobenius()`

Finds the Frobenius norm of \$matrix. This is just an alias of sqrt(\$matrix->sumeach(sub {\$_[0]*\$_[0]}).

• `\$matrix->trace()`

Finds the trace of \$matrix, which is the sum of its diagonal elements. This just is an alias for \$matrix->diagpart->sum().

• `\$matrix->diagnose()`

Sets the special flags of \$matrix.

## BOOLEAN INFORMATIONAL METHODS

The following is_xxx() methods use the special flags if they are present. If more than one will be called for a given matrix, consider calling \$matrix->diagnose() to set the flags.

• `\$matrix->is_square()`

`\$matrix->is_quadratic()`

Returns the value of the comparison \$matrix->{rows}==\$matrix->{columns} That is, it returns 1 if the matrix is square, 0 otherwise.

• `\$matrix->is_diagonal()`

Returns 1 if \$matrix is a diagonal matrix, 0 otherwise. \$matrix need not be square.

• `\$matrix->is_lowertriangular()`

Returns 1 if \$matrix is a lower triangular matrix, 0 otherwise. \$matrix need not be square. Diagonal elements are allowed.

• `\$matrix->is_strictlowertriangular()`

Returns 1 if \$matrix is a lower triangular matrix, 0 otherwise. \$matrix need not be square. Diagonal elements are not allowed.

• `\$matrix->is_uppertriangular()`

Returns 1 if \$matrix is an upper triangular matrix, 0 otherwise. \$matrix need not be square. Diagonal elements are allowed.

• `\$matrix->is_strictuppertriangular()`

Returns 1 if \$matrix is an upper triangular matrix, 0 otherwise. \$matrix need not be square. Diagonal elements are not allowed.

• `\$matrix->is_positive()`

Returns 1 if all elements of \$matrix are positive, 0 otherwise.

• `\$matrix->is_negative()`

Returns 1 if all elements of \$matrix are negative, 0 otherwise.

• `\$matrix->is_nonpositive()`

Returns 1 if all elements of \$matrix are nonpositive (i.e. <=0), 0 otherwise.

• `\$matrix->is_nonnegative()`

Returns 1 if all elements of \$matrix are nonnegative (i.e. >=0), 0 otherwise.

• `\$matrix->is_boolean()`

`\$matrix->is_pattern()`

Returns 1 if all elements of \$matrix are 0 or 1, 0 otherwise.

• `\$matrix->is_symmetric()`

Returns 1 if \$matrix is symmetric, i.e. \$matrix==\$matrix->transpose(), 0 otherwise. \$matrix is assumed to be square.

• `\$matrix->is_skewsymmetric()`

Returns 1 if \$matrix is skew-symmetric, i.e. \$matrix==-1*\$matrix->transpose(), 0 otherwise. \$matrix is assumed to be square.

## ARITHMETIC METHODS

• `\$matrix1->add(\$matrix2)`

Finds and returns \$matrix1 + \$matrix2. Matrices must be of the same size.

• `\$matrix1->multiply(\$matrix2)`

Finds and returns \$matrix1 * \$matrix2. \$matrix1->{columns} must equal \$matrix2->{rows}.

• `\$matrix1->multiplyfree(\$matrix2)`

`\$matrix1->addfree(\$matrix2)`

As add() and multiply(), but with no limitations on the dimensions of the matrices.

• `\$matrix1->quickmultiply(\$matrix2)`

`\$matrix1->quickmultiplyfree(\$matrix2)`

As multiply() and multiplyfree(), but with a different algorithm that is faster for larger matrices.

• `\$matrix->multiplyscalar(\$scalar)`

Returns the product of \$matrix and \$scalar.

• `\$matrix1->kronecker(\$matrix2)`

Returns the direct product of \$matrix1 and \$matrix2. Every element of \$matrix1 is multiplied by the entirely of \$matrix2, and those matrices assembled into a big matrix and returned.

• `\$matrix->exponentiate(\$n)` `\$matrix->largeexponentiate(\$n)`

Finds and returns \$matrix raised to the \$n th power. \$n must be nonnegative and integral, and \$matrix must be square.

For large values of \$n, largeexponentiate() is faster.

• `\$matrix->terminvert()`

Returns the matrix whose elements are the multiplicative inverses of \$matrix. If \$matrix is square diagonal, `\$matrix->terminvert()` is the multiplicative inverse of \$matrix.

• `\$matrix->transpose()`

Returns the transpose of \$matrix.

• `\$matrix->each(\$coderef)`

Applies a subroutine referred to by \$coderef to each element of \$matrix. &\$coderef should take three arguments (\$value, \$row, \$column).

• `\$matrix->symmetrify()`

Returns a matrix obtained by reflecting \$matrix about its main diagonal. \$matrix does not need to be square.

Exits with an error if \$matrix(i,j) and \$matrix(j,i) both exist and are not identical.

• `\$matrix->skewsymmetrify()`

As symmetrify(), except that the reflected terms are made negative.

Exits with an error if \$matrix(i,j) and \$matrix(j,i) both exist and are not negatives of each other.

• `\$matrix1->matrixand(\$matrix2)`

Finds and returns the matrix whose (\$i,\$j) element is \$matrix1(\$i,\$j) && \$matrix2(\$i,\$j). The returned matrix is a pattern matrix.

• `\$matrix1->matrixor(\$matrix2)`

Finds and returns the matrix whose (\$i,\$j) element is \$matrix1(\$i,\$j) || \$matrix2(\$i,\$j). The returned matrix is a pattern matrix.

## PROBABILISTIC METHODS

• `\$matrix->normalize()`

Returns a matrix \$matrix2 scaled so that \$matrix2->sum()==1.

Exits on error if \$matrix is not nonnegative.

• `\$matrix->normalizerows()`

As \$matrix->normalize(), except that each row is scaled to have a sum of 1.

Exits on error if \$matrix is not nonnegative.

• `\$matrix->normalizecolumns()`

As \$matrix->normalizerows(), except with columns. Mathematically equivalent to \$matrix->transpose()->normalizerows()->transpose().

Exits on error if \$matrix is not nonnegative.

• `\$matrix->discretepdf()`

Assuming that \$matrix->sum()==1 and \$matrix is non-negative, chooses a random element from \$matrix based on the assumption that the entry at (\$i,\$j) is the probability of choosing (\$i,\$j).

## SOLUTION OF SYSTEMS METHODS

• `\$matrix->jacobi(\$constant,\$guess, \$tol, \$steps)`

Uses Jacobi iteration to find and return the solution to the system of equations \$matrix * x = \$constant, with initial guess \$constant, tolerance \$tol, and maximum iterations \$steps.

If \$steps is undefined, the default value of 100 is used.

Care should be taken to ensure that \$matrix is such that the iteration actually converges.

• `\$matrix->gaussseidel(\$constant,\$guess, \$tol, \$steps)`

Uses Gauss-Seidel iteration to find and return the solution to the system of equations \$matrix * x = \$constant, with initial guess \$constant, tolerance \$tol, and maximum iterations \$steps. This is equivalent to \$matrix->SOR with relaxation parameter 1.

Care should be taken to ensure that \$matrix is such that the iteration actually converges.

• `\$matrix->SOR(\$constant,\$guess, \$relax, \$tol, \$steps)`

Uses Successive Over-Relaxation to find and return the solution to the system of equations \$matrix * x = \$constant, with initial guess \$constant, relaxation parameter \$relax, tolerance \$tol, and maximum iterations \$steps.

Care should be taken to ensure that \$matrix is such that the iteration actually converges.

## SORTING METHODS

• `\$matrix->sortbycolumn()`

`\$matrix->sortbyrow()`

Returns an array containing the keys of \$matrix, sorted primarily by either column or row (and secondarily by row or column).

`\$matrix->sortbydiagonal()`

Returns an array containing the keys of \$matrix, sorted primarily by their distance from elements on the main diagonal, lower diagonals first. The row of the key is a secondary criterion.

`\$matrix->sortbyvalue()`

Returns an array containing the keys of \$matrix, sorted primarily by the value of the element indexed by the key. Row and column are secondary and tertiaty criteria.

These methods are suitable for using inside a loop. See also \$matrix->elements() if the order of the elements is not important.

Example:

A 3x3 matrix will be sorted in the following manners:

``` sortbycolumn()
1 4 7
2 5 8
3 6 9

sortbyrow()
1 2 3
4 5 6
7 8 9

sortbydiagonal()
4 7 9
2 5 8
1 3 6```

## IN-LINE METHODS

The following are as described above, except that they modify their calling object instead of a copy thereof.

`_each`

`_mask`

`_insert`

`_cull`

`_threshold`

`_sizetofit`

`_negate`

`_multiplyscalar`

`_terminvert`

`_symmetrify`

`_skewsymmetrify`

`_normalize()`

`_normalizerows()`

`_normalizecolumns()`

`_diagpart()`

`_nondiagpart()`

`_upperpart()`

`_nonupperpart()`

`_lowerpart()`

`_nonlowerpart()`

`_positivepart()`

`_negativepart()`

`_symmetricpart()`

`_skewsymmetricpart()`

In addition to these, the following methods modify the calling object.

`name()`

`assign()`

`assignspecial()`

`assignkey()`

`row()` (if called with the optional second argument)

`column()` (if called with the optional second argument)

`diagnose()`

`delete()`

`deletekey()`

`clearspecials()`

`+` add()

`-` subtract()

`*` quickmultiply()

`x` kronecker()

`**` exponential()

`~` transpose()

`&` matrixand()

`|` matrixor()

`""` print()

Unary `-` negate()

## SPECIAL MATRIX FLAGS

Certain information is stored about the matrix, and is updated when necessary. See also `\$matrix->diagnose()`. All flags can also be `undef`. These should not be accessed directly--use the boolean is_whatever() methods instead.

structure

The symmetry of the matrix. Can be `symmetric` or `skewsymmetric`.

shape

The shape of the matrix. Can be `diagonal`, `lower` (triangular), `upper`, `strictlower` or `strictupper`.

sign

The sign of all the elemetns of a matrix. Can be `positive`, `negative`, `nonpostive`, `nonnegative`, or `zero`.

pattern

Indicates whether the matrix should be treated as a pattern. Is non-zero if so.

bandwidth

Calculates the bandwidth of \$matrix, i.e. the maximum value of abs(\$i-\$j) for all elements at (\$i,\$j) in \$matrix.

field

The underlying field of the elements of the matrix. Currently, can only be `real`.

## INTERFACING

The user should not attempt to access the pieces of a Math::MatrixReal object directly, but instead use the routines provided. Certain methods, such as the sorters, return keys to the elements of the matrix, and these should be fed into `splitkey()` to obtain row-column indices.

None.

# EXPORT_OK

&splitkey() and &makekey().

# BUGS

Horribly, hideously inefficient.

No checks are made to be sure that values are of a proper type, or even that indices are integers. It is entirely possible to assign() a value that is, say, another Math::MatrixSparse. However, because of the lack of these checks, indices can start at zero, or even negative values, if an algorithm calls for it.

Harwell Boeing support is not robust, and output is not implemented.

Complex matrices in Harwell Boeing and Matrix Market are not supported. In Matrix Market format, only the real part is read--the imaginary part is discarded.

Many methods do not modify their calling object, but instead return a new object. This can be inefficent and a waste of resources, especially when it will be assigned to the new object anyway. Use the analogous methods listed in IN-LINE METHODS instead if this is an issue.

# AUTHOR

Jacob C. Kesinger, <kesinger@math.ttu.edu>