Sparse Matrices
***************

The Creation and Manipulation of Sparse Matrices
================================================

The size of mathematical problems that can be treated at any particular time is generally limited by the available computing resources.  Both, the speed of the computer and its available memory place limitation on the problem size.

There are many classes of mathematical problems which give rise to matrices, where a large number of the elements are zero.  In this case it makes sense to have a special matrix type to handle this class of problems where only the non-zero elements of the matrix are stored.  Not only does this reduce the amount of memory to store the matrix, but it also means that operations on this type of matrix can take advantage of the a-priori knowledge of the positions of the non-zero elements to accelerate their calculations.

A matrix type that stores only the non-zero elements is generally called sparse.  It is the purpose of this document to discuss the basics of the storage and creation of sparse matrices and the fundamental operations on them.

Storage of Sparse Matrices
--------------------------

It is not strictly speaking necessary for the user to understand how sparse matrices are stored.  However, such an understanding will help to get an understanding of the size of sparse matrices.  Understanding the storage technique is also necessary for those users wishing to create their own oct-files.

There are many different means of storing sparse matrix data.  What all of the methods have in common is that they attempt to reduce the complexity and storage given a-priori knowledge of the particular class of problems that will be solved.  A good summary of the available techniques for storing sparse matrix is given by Saad  (Youcef Saad "SPARSKIT: A basic toolkit for sparse matrix computation", 1994, http://www-users.cs.umn.edu/~saad/software/SPARSKIT/paper.ps). With full matrices, knowledge of the point of an element of the matrix within the matrix is implied by its position in the computers memory. However, this is not the case for sparse matrices, and so the positions of the non-zero elements of the matrix must equally be stored.

An obvious way to do this is by storing the elements of the matrix as triplets, with two elements being their position in the array (rows and column) and the third being the data itself.  This is conceptually easy to grasp, but requires more storage than is strictly needed.

The storage technique used within Octave is the compressed column format.  In this format the position of each element in a row and the data are stored as previously.  However, if we assume that all elements in the same column are stored adjacent in the computers memory, then we only need to store information on the number of non-zero elements in each column, rather than their positions.  Thus assuming that the matrix has more non-zero elements than there are columns in the matrix, we win in terms of the amount of memory used.

In fact, the column index contains one more element than the number of columns, with the first element always being zero.  The advantage of this is a simplification in the code, in that there is no special case for the first or last columns.  A short example, demonstrating this in C is.

       for (j = 0; j < nc; j++)
         for (i = cidx (j); i < cidx(j+1); i++)
            printf ("non-zero element (%i,%i) is %d\n",
     	   ridx(i), j, data(i));

A clear understanding might be had by considering an example of how the above applies to an example matrix.  Consider the matrix

         1   2   0  0
         0   0   0  3
         0   0   0  4

The non-zero elements of this matrix are

        (1, 1)  => 1
        (1, 2)  => 2
        (2, 4)  => 3
        (3, 4)  => 4

This will be stored as three vectors CIDX, RIDX and DATA, representing the column indexing, row indexing and data respectively.  The contents of these three vectors for the above matrix will be

       CIDX = [0, 1, 2, 2, 4]
       RIDX = [0, 0, 1, 2]
       DATA = [1, 2, 3, 4]

Note that this is the representation of these elements with the first row and column assumed to start at zero, while in Octave itself the row and column indexing starts at one.  Thus the number of elements in the I-th column is given by CIDX (I + 1) - CIDX (I).

Although Octave uses a compressed column format, it should be noted that compressed row formats are equally possible.  However, in the context of mixed operations between mixed sparse and dense matrices, it makes sense that the elements of the sparse matrices are in the same order as the dense matrices.  Octave stores dense matrices in column major ordering, and so sparse matrices are equally stored in this manner.

A further constraint on the sparse matrix storage used by Octave is that all elements in the rows are stored in increasing order of their row index, which makes certain operations faster.  However, it imposes the need to sort the elements on the creation of sparse matrices.  Having disordered elements is potentially an advantage in that it makes operations such as concatenating two sparse matrices together easier and faster, however it adds complexity and speed problems elsewhere.

Creating Sparse Matrices
------------------------

There are several means to create sparse matrix.

Returned from a function
     There are many functions that directly return sparse matrices.  These include speye, sprand, diag, etc.

Constructed from matrices or vectors
     The function sparse allows a sparse matrix to be constructed from three vectors representing the row, column and data.  Alternatively, the function spconvert uses a three column matrix format to allow easy importation of data from elsewhere.

Created and then filled
     The function sparse or spalloc can be used to create an empty matrix that is then filled by the user

From a user binary program
     The user can directly create the sparse matrix within an oct-file.

There are several basic functions to return specific sparse matrices.  For example the sparse identity matrix, is a matrix that is often needed.  It therefore has its own function to create it as speye (N) or speye (R, C), which creates an N-by-N or R-by-C sparse identity matrix.

Another typical sparse matrix that is often needed is a random distribution of random elements.  The functions sprand and sprandn perform this for uniform and normal random distributions of elements.  They have exactly the same calling convention, where sprand (R, C, D), creates an R-by-C sparse matrix with a density of filled elements of D.

Other functions of interest that directly create sparse matrices, are diag or its generalization spdiags, that can take the definition of the diagonals of the matrix and create the sparse matrix that corresponds to this.  For example

     s = diag (sparse(randn(1,n)), -1);

creates a sparse (N+1)-by-(N+1) sparse matrix with a single diagonal defined.

 -- Function File: [B, C] = spdiags (A)
 -- Function File: B = spdiags (A, C)
 -- Function File: B = spdiags (V, C, A)
 -- Function File: B = spdiags (V, C, M, N)
     A generalization of the function diag.  Called with a single input argument, the non-zero diagonals C of A are extracted. With two arguments the diagonals to extract are given by the vector C.

     The other two forms of spdiags modify the input matrix by replacing the diagonals.  They use the columns of V to replace the columns represented by the vector C.  If the sparse matrix A is defined then the diagonals of this matrix are replaced. Otherwise a matrix of M by N is created with the diagonals given by V.

     Negative values of C represent diagonals below the main diagonal, and positive values of C diagonals above the main diagonal.

     For example

          spdiags (reshape (1:12, 4, 3), [-1 0 1], 5, 4)
          =>    5 10  0  0
                1  6 11  0
                0  2  7 12
                0  0  3  8
                0  0  0  4

 -- Function File: Y = speye (M)
 -- Function File: Y = speye (M, N)
 -- Function File: Y = speye (SZ)
     Returns a sparse identity matrix.  This is significantly more efficient than sparse (eye (M)) as the full matrix is not constructed.

     Called with a single argument a square matrix of size M by M is created.  Otherwise a matrix of M by N is created.  If called with a single vector argument, this argument is taken to be the size of the matrix to create.

 -- Function File: Y = spfun (F,X)
     Compute f(X) for the non-zero values of X. This results in a sparse matrix with the same structure as X.  The function F can be passed as a string, a function handle or an inline function.

 -- Mapping Function: spmax (X, Y, DIM)
 -- Mapping Function: [W, IW] = spmax (X)
     This function has been deprecated.  Use max instead.

 -- Mapping Function: spmin (X, Y, DIM)
 -- Mapping Function: [W, IW] = spmin (X)
     This function has been deprecated.  Use min instead.

 -- Function File: Y = spones (X)
     Replace the non-zero entries of X with ones.  This creates a sparse matrix with the same structure as X.

 -- Function File: sprand (M, N, D)
 -- Function File: sprand (S)
     Generate a random sparse matrix.  The size of the matrix will be M by N, with a density of values given by D. D should be between 0 and 1. Values will be uniformly distributed between 0 and 1.

     Note: sometimes the actual density may be a bit smaller than D. This is unlikely to happen for large really sparse matrices.

     If called with a single matrix argument, a random sparse matrix is generated wherever the matrix S is non-zero. sprandn

 -- Function File: sprandn (M, N, D)
 -- Function File: sprandn (S)
     Generate a random sparse matrix.  The size of the matrix will be M by N, with a density of values given by D. D should be between 0 and 1. Values will be normally distributed with mean of zero and variance 1.

     Note: sometimes the actual density may be a bit smaller than D. This is unlikely to happen for large really sparse matrices.

     If called with a single matrix argument, a random sparse matrix is generated wherever the matrix S is non-zero. sprand

 -- Function File: sprandsym (N, D)
 -- Function File: sprandsym (S)
     Generate a symmetric random sparse matrix.  The size of the matrix will be N by N, with a density of values given by D. D should be between 0 and 1. Values will be normally distributed with mean of zero and variance 1.

     Note: sometimes the actual density may be a bit smaller than D. This is unlikely to happen for large really sparse matrices.

     If called with a single matrix argument, a random sparse matrix is generated wherever the matrix S is non-zero in its lower triangular part. sprand, sprandn

The recommended way for the user to create a sparse matrix, is to create two vectors containing the row and column index of the data and a third vector of the same size containing the data to be stored.  For example

       ri = ci = d = [];
       for j = 1:c
         ri = [ri; randperm(r)(1:n)'];
         ci = [ci; j*ones(n,1)];
         d = [d; rand(n,1)];
       endfor
       s = sparse (ri, ci, d, r, c);

creates an R-by-C sparse matrix with a random distribution of N (<R) elements per column.  The elements of the vectors do not need to be sorted in any particular order as Octave will sort them prior to storing the data.  However, pre-sorting the data will make the creation of the sparse matrix faster.

The function spconvert takes a three or four column real matrix. The first two columns represent the row and column index respectively and the third and four columns, the real and imaginary parts of the sparse matrix.  The matrix can contain zero elements and the elements can be sorted in any order.  Adding zero elements is a convenient way to define the size of the sparse matrix.  For example

     s = spconvert ([1 2 3 4; 1 3 4 4; 1 2 3 0]')
     => Compressed Column Sparse (rows=4, cols=4, nnz=3)
           (1 , 1) -> 1
           (2 , 3) -> 2
           (3 , 4) -> 3

An example of creating and filling a matrix might be

     k = 5;
     nz = r * k;
     s = spalloc (r, c, nz)
     for j = 1:c
       idx = randperm (r);
       s (:, j) = [zeros(r - k, 1); ...
             rand(k, 1)] (idx);
     endfor

It should be noted, that due to the way that the Octave assignment functions are written that the assignment will reallocate the memory used by the sparse matrix at each iteration of the above loop. Therefore the spalloc function ignores the NZ argument and does not preassign the memory for the matrix.  Therefore, it is vitally important that code using to above structure should be vectorized as much as possible to minimize the number of assignments and reduce the number of memory allocations.

 -- Loadable Function: FM = full (SM)
     returns a full storage matrix from a sparse, diagonal, permutation matrix or a range. sparse

 -- Function File: S = spalloc (R, C, NZ)
     Returns an empty sparse matrix of size R-by-C.  As Octave resizes sparse matrices at the first opportunity, so that no additional space is needed, the argument NZ is ignored.  This function is provided only for compatibility reasons.

     It should be noted that this means that code like

          k = 5;
          nz = r * k;
          s = spalloc (r, c, nz)
          for j = 1:c
            idx = randperm (r);
            s (:, j) = [zeros(r - k, 1); rand(k, 1)] (idx);
          endfor

     will reallocate memory at each step.  It is therefore vitally important that code like this is vectorized as much as possible. sparse, nzmax

 -- Loadable Function: S = sparse (A)
 -- Loadable Function: S = sparse (I, J, SV, M, N, NZMAX)
 -- Loadable Function: S = sparse (I, J, SV)
 -- Loadable Function: S = sparse (I, J, S, M, N, "unique")
 -- Loadable Function: S = sparse (M, N)
     Create a sparse matrix from the full matrix or row, column, value triplets. If A is a full matrix, convert it to a sparse matrix representation, removing all zero values in the process.

     Given the integer index vectors I and J, a 1-by-nnz vector of real of complex values SV, overall dimensions M and N of the sparse matrix.  The argument nzmax is ignored but accepted for compatibility with MATLAB.  If M or N are not specified their values are derived from the maximum index in the vectors I and J as given by M = max (I), N = max (J).

     Note: if multiple values are specified with the same I, J indices, the corresponding values in S will be added.

     The following are all equivalent:

          s = sparse (i, j, s, m, n)
          s = sparse (i, j, s, m, n, "summation")
          s = sparse (i, j, s, m, n, "sum")

     Given the option "unique". if more than two values are specified for the same I, J indices, the last specified value will be used.

     sparse(M, N) is equivalent to sparse ([], [], [], M, N, 0)

     If any of SV, I or J are scalars, they are expanded to have a common size. full

 -- Function File: X = spconvert (M)
     This function converts for a simple sparse matrix format easily produced by other programs into Octave's internal sparse format.  The input X is either a 3 or 4 column real matrix, containing the row, column, real and imaginary parts of the elements of the sparse matrix.  An element with a zero real and imaginary part can be used to force a particular matrix size.

The above problem of memory reallocation can be avoided in oct-files.  However, the construction of a sparse matrix from an oct-file is more complex than can be discussed here, and you are referred to chapter Dynamically Linked Functions, to have a full description of the techniques involved.

Finding out Information about Sparse Matrices
---------------------------------------------

There are a number of functions that allow information concerning sparse matrices to be obtained.  The most basic of these is issparse that identifies whether a particular Octave object is in fact a sparse matrix.

Another very basic function is nnz that returns the number of non-zero entries there are in a sparse matrix, while the function nzmax returns the amount of storage allocated to the sparse matrix.  Note that Octave tends to crop unused memory at the first opportunity for sparse objects.  There are some cases of user created sparse objects where the value returned by nzmax will not be the same as nnz, but in general they will give the same result.  The function spstats returns some basic statistics on the columns of a sparse matrix including the number of elements, the mean and the variance of each column.

 -- Loadable Function: issparse (EXPR)
     Return 1 if the value of the expression EXPR is a sparse matrix.

 -- Built-in Function: SCALAR = nnz (A)
     Returns the number of non zero elements in A. sparse

 -- Function File: nonzeros (S)
     Returns a vector of the non-zero values of the sparse matrix S.

 -- Built-in Function: SCALAR = nzmax (SM)
     Return the amount of storage allocated to the sparse matrix SM. Note that Octave tends to crop unused memory at the first opportunity for sparse objects.  There are some cases of user created sparse objects where the value returned by nzmax will not be the same as nnz, but in general they will give the same result. sparse, spalloc

 -- Function File: [COUNT, MEAN, VAR] = spstats (S)
 -- Function File: [COUNT, MEAN, VAR] = spstats (S, J)
     Return the stats for the non-zero elements of the sparse matrix S. COUNT is the number of non-zeros in each column, MEAN is the mean of the non-zeros in each column, and VAR is the variance of the non-zeros in each column.

     Called with two input arguments, if S is the data and J is the bin number for the data, compute the stats for each bin.  In this case, bins can contain data values of zero, whereas with spstats (S) the zeros may disappear.

When solving linear equations involving sparse matrices Octave determines the means to solve the equation based on the type of the matrix as discussed in Sparse Linear Algebra.  Octave probes the matrix type when the div (/) or ldiv (\) operator is first used with the matrix and then caches the type.  However the matrix_type function can be used to determine the type of the sparse matrix prior to use of the div or ldiv operators.  For example

     a = tril (sprandn(1024, 1024, 0.02), -1) ...
         + speye(1024);
     matrix_type (a);
     ans = Lower

show that Octave correctly determines the matrix type for lower triangular matrices.  matrix_type can also be used to force the type of a matrix to be a particular type.  For example

     a = matrix_type (tril (sprandn (1024, ...
        1024, 0.02), -1) + speye(1024), 'Lower');

This allows the cost of determining the matrix type to be avoided.  However, incorrectly defining the matrix type will result in incorrect results from solutions of linear equations, and so it is entirely the responsibility of the user to correctly identify the matrix type

There are several graphical means of finding out information about sparse matrices.  The first is the spy command, which displays the structure of the non-zero elements of the matrix.  See fig:spmatrix, for an example of the use of spy.  More advanced graphical information can be obtained with the treeplot, etreeplot and gplot commands.

Figure: Structure of simple sparse matrix.

One use of sparse matrices is in graph theory, where the interconnections between nodes are represented as an adjacency matrix.  That is, if the i-th node in a graph is connected to the j-th node.  Then the ij-th node (and in the case of undirected graphs the ji-th node) of the sparse adjacency matrix is non-zero.  If each node is then associated with a set of coordinates, then the gplot command can be used to graphically display the interconnections between nodes.

As a trivial example of the use of gplot, consider the example

     A = sparse([2,6,1,3,2,4,3,5,4,6,1,5],
         [1,1,2,2,3,3,4,4,5,5,6,6],1,6,6);
     xy = [0,4,8,6,4,2;5,0,5,7,5,7]';
     gplot(A,xy)

which creates an adjacency matrix A where node 1 is connected to nodes 2 and 6, node 2 with nodes 1 and 3, etc.  The coordinates of the nodes are given in the n-by-2 matrix xy.

The dependencies between the nodes of a Cholesky factorization can be calculated in linear time without explicitly needing to calculate the Cholesky factorization by the etree command.  This command returns the elimination tree of the matrix and can be displayed graphically by the command treeplot(etree(A)) if A is symmetric or treeplot(etree(A+A')) otherwise.

 -- Function File: spy (X)
 -- Function File: spy (..., MARKERSIZE)
 -- Function File: spy (..., LINE_SPEC)
     Plot the sparsity pattern of the sparse matrix X.  If the argument MARKERSIZE is given as an scalar value, it is used to determine the point size in the plot.  If the string LINE_SPEC is given it is passed to plot and determines the appearance of the plot. plot

 -- Loadable Function: P = etree (S)
 -- Loadable Function: P = etree (S, TYP)
 -- Loadable Function: [P, Q] = etree (S, TYP)

     Returns the elimination tree for the matrix S.  By default S is assumed to be symmetric and the symmetric elimination tree is returned.  The argument TYP controls whether a symmetric or column elimination tree is returned.  Valid values of TYP are 'sym' or 'col', for symmetric or column elimination tree respectively

     Called with a second argument, etree also returns the postorder permutations on the tree.

 -- Function File: etreeplot (TREE)
 -- Function File: etreeplot (TREE, NODE_STYLE, EDGE_STYLE)
     Plot the elimination tree of the matrix S or S+S'  if S in non-symmetric.  The optional parameters LINE_STYLE and EDGE_STYLE define the output style. treeplot, gplot

 -- Function File: gplot (A, XY)
 -- Function File: gplot (A, XY, LINE_STYLE)
 -- Function File: [X, Y] = gplot (A, XY)
     Plot a graph defined by A and XY in the graph theory sense.  A is the adjacency matrix of the array to be plotted and XY is an N-by-2 matrix containing the coordinates of the nodes of the graph.

     The optional parameter LINE_STYLE defines the output style for the plot.  Called with no output arguments the graph is plotted directly.  Otherwise, return the coordinates of the plot in X and Y. treeplot, etreeplot, spy

 -- Function File: treeplot (TREE)
 -- Function File: treeplot (TREE, LINE_STYLE, EDGE_STYLE)
     Produces a graph of tree or forest.  The first argument is vector of predecessors, optional parameters LINE_STYLE and EDGE_STYLE define the output style.  The complexity of the algorithm is O(n) in terms of is time and memory requirements. etreeplot, gplot

 -- Function File: treelayout (TREE)
 -- Function File: treelayout (TREE, PERMUTATION)
     treelayout lays out a tree or a forest.  The first argument TREE is a vector of predecessors, optional parameter PERMUTATION is an optional postorder permutation. The complexity of the algorithm is O(n) in terms of time and memory requirements. etreeplot, gplot, treeplot

Basic Operators and Functions on Sparse Matrices
------------------------------------------------

Sparse Functions
................

An important consideration in the use of the sparse functions of Octave is that many of the internal functions of Octave, such as diag, cannot accept sparse matrices as an input.  The sparse implementation in Octave therefore uses the dispatch function to overload the normal Octave functions with equivalent functions that work with sparse matrices.  However, at any time the sparse matrix specific version of the function can be used by explicitly calling its function name.

The table below lists all of the sparse functions of Octave.  Note that the names of the specific sparse forms of the functions are typically the same as the general versions with a sp prefix.  In the table below, and the rest of this article the specific sparse versions of the functions are used.

Generate sparse matrices:
     spalloc, spdiags, speye, sprand, sprandn, sprandsym

Sparse matrix conversion:
     full, sparse, spconvert

Manipulate sparse matrices
     issparse, nnz, nonzeros, nzmax, spfun, spones, spy

Graph Theory:
     etree, etreeplot, gplot, treeplot

Sparse matrix reordering:
     amd, ccolamd, colamd, colperm, csymamd, dmperm, symamd, randperm, symrcm

Linear algebra:
     condest, eigs, matrix_type, normest, sprank, spaugment, svds

Iterative techniques:
     luinc, pcg, pcr

Miscellaneous:
     spparms, symbfact, spstats

In addition all of the standard Octave mapper functions (i.e., basic math functions that take a single argument) such as abs, etc. can accept sparse matrices.  The reader is referred to the documentation supplied with these functions within Octave itself for further details.

The Return Types of Operators and Functions
...........................................

The two basic reasons to use sparse matrices are to reduce the memory usage and to not have to do calculations on zero elements.  The two are closely related in that the computation time on a sparse matrix operator or function is roughly linear with the number of non-zero elements.

Therefore, there is a certain density of non-zero elements of a matrix where it no longer makes sense to store it as a sparse matrix, but rather as a full matrix.  For this reason operators and functions that have a high probability of returning a full matrix will always return one.  For example adding a scalar constant to a sparse matrix will almost always make it a full matrix, and so the example

     speye(3) + 0
     =>   1  0  0
       0  1  0
       0  0  1

returns a full matrix as can be seen.

Additionally, if sparse_auto_mutate is true, all sparse functions test the amount of memory occupied by the sparse matrix to see if the amount of storage used is larger than the amount used by the full equivalent.  Therefore speye (2) * 1 will return a full matrix as the memory used is smaller for the full version than the sparse version.

As all of the mixed operators and functions between full and sparse matrices exist, in general this does not cause any problems.  However, one area where it does cause a problem is where a sparse matrix is promoted to a full matrix, where subsequent operations would resparsify the matrix.  Such cases are rare, but can be artificially created, for example (fliplr(speye(3)) + speye(3)) - speye(3) gives a full matrix when it should give a sparse one.  In general, where such cases occur, they impose only a small memory penalty.

There is however one known case where this behavior of Octave's sparse matrices will cause a problem.  That is in the handling of the diag function.  Whether diag returns a sparse or full matrix depending on the type of its input arguments.  So

      a = diag (sparse([1,2,3]), -1);

should return a sparse matrix.  To ensure this actually happens, the sparse function, and other functions based on it like speye, always returns a sparse matrix, even if the memory used will be larger than its full representation.

 -- Built-in Function: VAL = sparse_auto_mutate ()
 -- Built-in Function: OLD_VAL = sparse_auto_mutate (NEW_VAL)
     Query or set the internal variable that controls whether Octave will automatically mutate sparse matrices to real matrices to save memory. For example,

          s = speye(3);
          sparse_auto_mutate (false)
          s (:, 1) = 1;
          typeinfo (s)
          => sparse matrix
          sparse_auto_mutate (true)
          s (1, :) = 1;
          typeinfo (s)
          => matrix

Note that the sparse_auto_mutate option is incompatible with MATLAB, and so it is off by default.

Mathematical Considerations
...........................

The attempt has been made to make sparse matrices behave in exactly the same manner as there full counterparts.  However, there are certain differences and especially differences with other products sparse implementations.

Firstly, the "./" and ".^" operators must be used with care.  Consider what the examples

       s = speye (4);
       a1 = s .^ 2;
       a2 = s .^ s;
       a3 = s .^ -2;
       a4 = s ./ 2;
       a5 = 2 ./ s;
       a6 = s ./ s;

will give.  The first example of S raised to the power of 2 causes no problems.  However S raised element-wise to itself involves a large number of terms 0 .^ 0 which is 1. There S .^ S is a full matrix.

Likewise S .^ -2 involves terms like 0 .^ -2 which is infinity, and so S .^ -2 is equally a full matrix.

For the "./" operator S ./ 2 has no problems, but 2 ./ S involves a large number of infinity terms as well and is equally a full matrix.  The case of S ./ S involves terms like 0 ./ 0 which is a NaN and so this is equally a full matrix with the zero elements of S filled with NaN values.

The above behavior is consistent with full matrices, but is not consistent with sparse implementations in other products.

A particular problem of sparse matrices comes about due to the fact that as the zeros are not stored, the sign-bit of these zeros is equally not stored.  In certain cases the sign-bit of zero is important.  For example

      a = 0 ./ [-1, 1; 1, -1];
      b = 1 ./ a
      => -Inf            Inf
          Inf           -Inf
      c = 1 ./ sparse (a)
      =>  Inf            Inf
          Inf            Inf

To correct this behavior would mean that zero elements with a negative sign-bit would need to be stored in the matrix to ensure that their sign-bit was respected.  This is not done at this time, for reasons of efficiency, and so the user is warned that calculations where the sign-bit of zero is important must not be done using sparse matrices.

In general any function or operator used on a sparse matrix will result in a sparse matrix with the same or a larger number of non-zero elements than the original matrix.  This is particularly true for the important case of sparse matrix factorizations.  The usual way to address this is to reorder the matrix, such that its factorization is sparser than the factorization of the original matrix.  That is the factorization of L * U = P * S * Q has sparser terms L and U than the equivalent factorization L * U = S.

Several functions are available to reorder depending on the type of the matrix to be factorized.  If the matrix is symmetric positive-definite, then symamd or csymamd should be used.  Otherwise amd, colamd or ccolamd should be used.  For completeness the reordering functions colperm and randperm are also available.

See fig:simplematrix, for an example of the structure of a simple positive definite matrix.

Figure: Structure of simple sparse matrix.

The standard Cholesky factorization of this matrix can be obtained by the same command that would be used for a full matrix.  This can be visualized with the command r = chol(A); spy(r);. The original matrix had 43 non-zero terms, while this Cholesky factorization has 71, with only half of the symmetric matrix being stored.  This is a significant level of fill in, and although not an issue for such a small test case, can represents a large overhead in working with other sparse matrices.

The appropriate sparsity preserving permutation of the original matrix is given by symamd and the factorization using this reordering can be visualized using the command q = symamd(A); r = chol(A(q,q)); spy(r).  This gives 29 non-zero terms which is a significant improvement.

The Cholesky factorization itself can be used to determine the appropriate sparsity preserving reordering of the matrix during the factorization, In that case this might be obtained with three return arguments as r[r, p, q] = chol(A); spy(r).

In the case of an asymmetric matrix, the appropriate sparsity preserving permutation is colamd and the factorization using this reordering can be visualized using the command q = colamd(A); [l, u, p] = lu(A(:,q)); spy(l+u).

Finally, Octave implicitly reorders the matrix when using the div (/) and ldiv (\) operators, and so no the user does not need to explicitly reorder the matrix to maximize performance.

 -- Loadable Function: P = amd (S)
 -- Loadable Function: P = amd (S, OPTS)

     Returns the approximate minimum degree permutation of a matrix.  This permutation such that the Cholesky factorization of S (P, P) tends to be sparser than the Cholesky factorization of S itself.  amd is typically faster than symamd but serves a similar purpose.

     The optional parameter OPTS is a structure that controls the behavior of amd.  The fields of these structure are

     opts.dense
          Determines what amd considers to be a dense row or column of the input matrix.  Rows or columns with more than max(16, (dense * sqrt (N) entries, where N is the order of the matrix S, are ignored by amd during the calculation of the permutation The value of dense must be a positive scalar and its default value is 10.0

     opts.aggressive
          If this value is a non zero scalar, then amd performs aggressive absorption.  The default is not to perform aggressive absorption.

     The author of the code itself is Timothy A. Davis (davis@cise.ufl.edu), University of Florida (see http://www.cise.ufl.edu/research/sparse/amd). symamd, colamd

 -- Loadable Function: P = ccolamd (S)
 -- Loadable Function: P = ccolamd (S, KNOBS)
 -- Loadable Function: P = ccolamd (S, KNOBS, CMEMBER)
 -- Loadable Function: [P, STATS] = ccolamd (...)

     Constrained column approximate minimum degree permutation.  P = ccolamd (S) returns the column approximate minimum degree permutation vector for the sparse matrix S.  For a non-symmetric matrix S, S (:, P) tends to have sparser LU factors than S. chol (S (:, P)' * S (:, P)) also tends to be sparser than chol (S' * S).  P = ccolamd (S, 1) optimizes the ordering for lu (S (:, P)). The ordering is followed by a column elimination tree post-ordering.

     KNOBS is an optional one- to five-element input vector, with a default value of [0 10 10 1 0] if not present or empty.  Entries not present are set to their defaults.

     KNOBS(1)
          if nonzero, the ordering is optimized for lu (S (:, p)).  It will be a poor ordering for chol (S (:, P)' * S (:, P)).  This is the most important knob for ccolamd.

     KNOB(2)
          if S is m-by-n, rows with more than max (16, KNOBS (2) * sqrt (n)) entries are ignored.

     KNOB(3)
          columns with more than max (16, KNOBS (3) * sqrt (min (M, N))) entries are ignored and ordered last in the output permutation (subject to the cmember constraints).

     KNOB(4)
          if nonzero, aggressive absorption is performed.

     KNOB(5)
          if nonzero, statistics and knobs are printed.

     CMEMBER is an optional vector of length n.  It defines the constraints on the column ordering.  If CMEMBER (j) = C, then column J is in constraint set C (C must be in the range 1 to N).  In the output permutation P, all columns in set 1 appear first, followed by all columns in set 2, and so on.  CMEMBER = ones(1,n) if not present or empty.  ccolamd (S, [], 1 : N) returns 1 : N

     P = ccolamd (S) is about the same as P = colamd (S).  KNOBS and its default values differ.  colamd always does aggressive absorption, and it finds an ordering suitable for both lu (S (:, P)) and chol (S (:, P)' * S (:, P)); it cannot optimize its ordering for lu (S (:, P)) to the extent that ccolamd (S, 1) can.

     STATS is an optional 20-element output vector that provides data about the ordering and the validity of the input matrix S.  Ordering statistics are in STATS (1 : 3).  STATS (1) and STATS (2) are the number of dense or empty rows and columns ignored by CCOLAMD and STATS (3) is the number of garbage collections performed on the internal data structure used by CCOLAMD (roughly of size 2.2 * nnz (S) + 4 * M + 7 * N integers).

     STATS (4 : 7) provide information if CCOLAMD was able to continue.  The matrix is OK if STATS (4) is zero, or 1 if invalid.  STATS (5) is the rightmost column index that is unsorted or contains duplicate entries, or zero if no such column exists. STATS (6) is the last seen duplicate or out-of-order row index in the column index given by STATS (5), or zero if no such row index exists.  STATS (7) is the number of duplicate or out-of-order row indices.  STATS (8 : 20) is always zero in the current version of CCOLAMD (reserved for future use).

     The authors of the code itself are S. Larimore, T. Davis (Uni of Florida) and S. Rajamanickam in collaboration with J. Bilbert and E. Ng.  Supported by the National Science Foundation (DMS-9504974, DMS-9803599, CCR-0203270), and a grant from Sandia National Lab.  See http://www.cise.ufl.edu/research/sparse for ccolamd, csymamd, amd, colamd, symamd, and other related orderings. colamd, csymamd

 -- Loadable Function: P = colamd (S)
 -- Loadable Function: P = colamd (S, KNOBS)
 -- Loadable Function: [P, STATS] = colamd (S)
 -- Loadable Function: [P, STATS] = colamd (S, KNOBS)

     Column approximate minimum degree permutation.  P = colamd (S) returns the column approximate minimum degree permutation vector for the sparse matrix S.  For a non-symmetric matrix S, S (:,P) tends to have sparser LU factors than S. The Cholesky factorization of S (:,P)' * S (:,P) also tends to be sparser than that of S' * S.

     KNOBS is an optional one- to three-element input vector.  If S is m-by-n, then rows with more than max(16,KNOBS(1)*sqrt(n)) entries are ignored.  Columns with more than max(16,knobs(2)*sqrt(min(m,n))) entries are removed prior to ordering, and ordered last in the output permutation P.  Only completely dense rows or columns are removed if KNOBS (1) and KNOBS (2) are < 0, respectively. If KNOBS (3) is nonzero, STATS and KNOBS are printed.  The default is KNOBS = [10 10 0].  Note that KNOBS differs from earlier versions of colamd

     STATS is an optional 20-element output vector that provides data about the ordering and the validity of the input matrix S.  Ordering statistics are in STATS (1:3).  STATS (1) and STATS (2) are the number of dense or empty rows and columns ignored by COLAMD and STATS (3) is the number of garbage collections performed on the internal data structure used by COLAMD (roughly of size 2.2 * nnz(S) + 4 * M + 7 * N integers).

     Octave built-in functions are intended to generate valid sparse matrices, with no duplicate entries, with ascending row indices of the nonzeros in each column, with a non-negative number of entries in each column (!) and so on.  If a matrix is invalid, then COLAMD may or may not be able to continue.  If there are duplicate entries (a row index appears two or more times in the same column) or if the row indices in a column are out of order, then COLAMD can correct these errors by ignoring the duplicate entries and sorting each column of its internal copy of the matrix S (the input matrix S is not repaired, however).  If a matrix is invalid in other ways then COLAMD cannot continue, an error message is printed, and no output arguments (P or STATS) are returned. COLAMD is thus a simple way to check a sparse matrix to see if it's valid.

     STATS (4:7) provide information if COLAMD was able to continue.  The matrix is OK if STATS (4) is zero, or 1 if invalid.  STATS (5) is the rightmost column index that is unsorted or contains duplicate entries, or zero if no such column exists. STATS (6) is the last seen duplicate or out-of-order row index in the column index given by STATS (5), or zero if no such row index exists.  STATS (7) is the number of duplicate or out-of-order row indices.  STATS (8:20) is always zero in the current version of COLAMD (reserved for future use).

     The ordering is followed by a column elimination tree post-ordering.

     The authors of the code itself are Stefan I. Larimore and Timothy A. Davis (davis@cise.ufl.edu), University of Florida.  The algorithm was developed in collaboration with John Gilbert, Xerox PARC, and Esmond Ng, Oak Ridge National Laboratory.  (see http://www.cise.ufl.edu/research/sparse/colamd) colperm, symamd

 -- Function File: P = colperm (S)
     Returns the column permutations such that the columns of S (:, P) are ordered in terms of increase number of non-zero elements.  If S is symmetric, then P is chosen such that S (P, P) orders the rows and columns with increasing number of non zeros elements.

 -- Loadable Function: P = csymamd (S)
 -- Loadable Function: P = csymamd (S, KNOBS)
 -- Loadable Function: P = csymamd (S, KNOBS, CMEMBER)
 -- Loadable Function: [P, STATS] = csymamd (...)

     For a symmetric positive definite matrix S, returns the permutation vector P such that S(P,P) tends to have a sparser Cholesky factor than S.  Sometimes csymamd works well for symmetric indefinite matrices too.  The matrix S is assumed to be symmetric; only the strictly lower triangular part is referenced. S must be square.  The ordering is followed by an elimination tree post-ordering.

     KNOBS is an optional one- to three-element input vector, with a default value of [10 1 0] if present or empty.  Entries not present are set to their defaults.

     KNOBS(1)
          If S is n-by-n, then rows and columns with more than max(16,KNOBS(1)*sqrt(n)) entries are ignored, and ordered last in the output permutation (subject to the cmember constraints).

     KNOBS(2)
          If nonzero, aggressive absorption is performed.

     KNOBS(3)
          If nonzero, statistics and knobs are printed.

     CMEMBER is an optional vector of length n. It defines the constraints on the ordering.  If CMEMBER(j) = S, then row/column j is in constraint set C (C must be in the range 1 to n).  In the output permutation P, rows/columns in set 1 appear first, followed by all rows/columns in set 2, and so on.  CMEMBER = ones(1,n) if not present or empty.  csymamd(S,[],1:n) returns 1:n.

     P = csymamd(S) is about the same as P = symamd(S).  KNOBS and its default values differ.

     STATS (4:7) provide information if CCOLAMD was able to continue.  The matrix is OK if STATS (4) is zero, or 1 if invalid.  STATS (5) is the rightmost column index that is unsorted or contains duplicate entries, or zero if no such column exists. STATS (6) is the last seen duplicate or out-of-order row index in the column index given by STATS (5), or zero if no such row index exists.  STATS (7) is the number of duplicate or out-of-order row indices.  STATS (8:20) is always zero in the current version of CCOLAMD (reserved for future use).

     The authors of the code itself are S. Larimore, T. Davis (Uni of Florida) and S. Rajamanickam in collaboration with J. Bilbert and E. Ng.  Supported by the National Science Foundation (DMS-9504974, DMS-9803599, CCR-0203270), and a grant from Sandia National Lab.  See http://www.cise.ufl.edu/research/sparse for ccolamd, csymamd, amd, colamd, symamd, and other related orderings. symamd, ccolamd

 -- Loadable Function: P = dmperm (S)
 -- Loadable Function: [P, Q, R, S] = dmperm (S)

     Perform a Dulmage-Mendelsohn permutation on the sparse matrix S. With a single output argument dmperm performs the row permutations P such that S (P,:) has no zero elements on the diagonal.

     Called with two or more output arguments, returns the row and column permutations, such that S (P, Q) is in block triangular form.  The values of R and S define the boundaries of the blocks.  If S is square then R == S.

     The method used is described in: A. Pothen & C.-J. Fan. Computing the block triangular form of a sparse matrix. ACM Trans. Math. Software, 16(4):303-324, 1990. colamd, ccolamd

 -- Loadable Function: P = symamd (S)
 -- Loadable Function: P = symamd (S, KNOBS)
 -- Loadable Function: [P, STATS] = symamd (S)
 -- Loadable Function: [P, STATS] = symamd (S, KNOBS)

     For a symmetric positive definite matrix S, returns the permutation vector p such that S (P, P) tends to have a sparser Cholesky factor than S.  Sometimes SYMAMD works well for symmetric indefinite matrices too.  The matrix S is assumed to be symmetric; only the strictly lower triangular part is referenced.  S must be square.

     KNOBS is an optional one- to two-element input vector.  If S is n-by-n, then rows and columns with more than max(16,KNOBS(1)*sqrt(n)) entries are removed prior to ordering, and ordered last in the output permutation P.  No rows/columns are removed if KNOBS(1) < 0.  If KNOBS (2) is nonzero, stats and KNOBS are printed.  The default is KNOBS = [10 0].  Note that KNOBS differs from earlier versions of symamd.

     STATS is an optional 20-element output vector that provides data about the ordering and the validity of the input matrix S.  Ordering statistics are in STATS (1:3).  STATS (1) = STATS (2) is the number of dense or empty rows and columns ignored by SYMAMD and STATS (3) is the number of garbage collections performed on the internal data structure used by SYMAMD (roughly of size 8.4 * nnz (tril (S, -1)) + 9 * N integers).

     Octave built-in functions are intended to generate valid sparse matrices, with no duplicate entries, with ascending row indices of the nonzeros in each column, with a non-negative number of entries in each column (!) and so on.  If a matrix is invalid, then SYMAMD may or may not be able to continue.  If there are duplicate entries (a row index appears two or more times in the same column) or if the row indices in a column are out of order, then SYMAMD can correct these errors by ignoring the duplicate entries and sorting each column of its internal copy of the matrix S (the input matrix S is not repaired, however).  If a matrix is invalid in other ways then SYMAMD cannot continue, an error message is printed, and no output arguments (P or STATS) are returned.  SYMAMD is thus a simple way to check a sparse matrix to see if it's valid.

     STATS (4:7) provide information if SYMAMD was able to continue.  The matrix is OK if STATS (4) is zero, or 1 if invalid.  STATS (5) is the rightmost column index that is unsorted or contains duplicate entries, or zero if no such column exists.  STATS (6) is the last seen duplicate or out-of-order row index in the column index given by STATS (5), or zero if no such row index exists.  STATS (7) is the number of duplicate or out-of-order row indices.  STATS (8:20) is always zero in the current version of SYMAMD (reserved for future use).

     The ordering is followed by a column elimination tree post-ordering.

     The authors of the code itself are Stefan I. Larimore and Timothy A. Davis (davis@cise.ufl.edu), University of Florida.  The algorithm was developed in collaboration with John Gilbert, Xerox PARC, and Esmond Ng, Oak Ridge National Laboratory.  (see http://www.cise.ufl.edu/research/sparse/colamd) colperm, colamd

 -- Loadable Function: P = symrcm (S)
     Symmetric reverse Cuthill-McKee permutation of S. Return a permutation vector P such that S (P, P) tends to have its diagonal elements closer to the diagonal than S.  This is a good preordering for LU or Cholesky factorization of matrices that come from 'long, skinny' problems.  It works for both symmetric and asymmetric S.

     The algorithm represents a heuristic approach to the NP-complete bandwidth minimization problem.  The implementation is based in the descriptions found in

     E. Cuthill, J. McKee: Reducing the Bandwidth of Sparse Symmetric Matrices. Proceedings of the 24th ACM National Conference, 157--172 1969, Brandon Press, New Jersey.

     Alan George, Joseph W. H. Liu: Computer Solution of Large Sparse Positive Definite Systems, Prentice Hall Series in Computational Mathematics, ISBN 0-13-165274-5, 1981.

     colperm, colamd, symamd

Linear Algebra on Sparse Matrices
=================================

Octave includes a polymorphic solver for sparse matrices, where the exact solver used to factorize the matrix, depends on the properties of the sparse matrix itself.  Generally, the cost of determining the matrix type is small relative to the cost of factorizing the matrix itself, but in any case the matrix type is cached once it is calculated, so that it is not re-determined each time it is used in a linear equation.

The selection tree for how the linear equation is solve is

 1. If the matrix is diagonal, solve directly and goto 8

 2. If the matrix is a permuted diagonal, solve directly taking into account the permutations.  Goto 8

 3. If the matrix is square, banded and if the band density is less than that given by spparms ("bandden") continue, else goto 4.

      a. If the matrix is tridiagonal and the right-hand side is not sparse continue, else goto 3b.

           1. If the matrix is hermitian, with a positive real diagonal, attempt Cholesky factorization using LAPACK xPTSV.

           2. If the above failed or the matrix is not hermitian with a positive real diagonal use Gaussian elimination with pivoting using LAPACK xGTSV, and goto 8.

      b. If the matrix is hermitian with a positive real diagonal, attempt Cholesky factorization using LAPACK xPBTRF.

      c. if the above failed or the matrix is not hermitian with a positive real diagonal use Gaussian elimination with pivoting using LAPACK xGBTRF, and goto 8.

 4. If the matrix is upper or lower triangular perform a sparse forward or backward substitution, and goto 8

 5. If the matrix is a upper triangular matrix with column permutations or lower triangular matrix with row permutations, perform a sparse forward or backward substitution, and goto 8

 6. If the matrix is square, hermitian with a real positive diagonal, attempt sparse Cholesky factorization using CHOLMOD.

 7. If the sparse Cholesky factorization failed or the matrix is not hermitian with a real positive diagonal, and the matrix is square, factorize using UMFPACK.

 8. If the matrix is not square, or any of the previous solvers flags a singular or near singular matrix, find a minimum norm solution using CXSPARSE (The CHOLMOD, UMFPACK and CXSPARSE packages were written by Tim Davis and are available at http://www.cise.ufl.edu/research/sparse/).

The band density is defined as the number of non-zero values in the matrix divided by the number of non-zero values in the matrix.  The banded matrix solvers can be entirely disabled by using spparms to set bandden to 1 (i.e., spparms ("bandden", 1)).

The QR solver factorizes the problem with a Dulmage-Mendelsohn, to separate the problem into blocks that can be treated as over-determined, multiple well determined blocks, and a final over-determined block.  For matrices with blocks of strongly connected nodes this is a big win as LU decomposition can be used for many blocks.  It also significantly improves the chance of finding a solution to over-determined problems rather than just returning a vector of NaN's.

All of the solvers above, can calculate an estimate of the condition number.  This can be used to detect numerical stability problems in the solution and force a minimum norm solution to be used.  However, for narrow banded, triangular or diagonal matrices, the cost of calculating the condition number is significant, and can in fact exceed the cost of factoring the matrix.  Therefore the condition number is not calculated in these cases, and Octave relies on simpler techniques to detect singular matrices or the underlying LAPACK code in the case of banded matrices.

The user can force the type of the matrix with the matrix_type function.  This overcomes the cost of discovering the type of the matrix. However, it should be noted that identifying the type of the matrix incorrectly will lead to unpredictable results, and so matrix_type should be used with care.

 -- Function File: [N, C] = normest (A, TOL)
     Estimate the 2-norm of the matrix A using a power series analysis.  This is typically used for large matrices, where the cost of calculating the norm (A) is prohibitive and an approximation to the 2-norm is acceptable.

     TOL is the tolerance to which the 2-norm is calculated.  By default TOL is 1e-6.  C returns the number of iterations needed for normest to converge.

 -- Function File: [EST, V, W, ITER] = onenormest (A, T)
 -- Function File: [EST, V, W, ITER] = onenormest (APPLY, APPLY_T, N, T)

     Apply Higham and Tisseur's randomized block 1-norm estimator to matrix A using T test vectors.  If T exceeds 5, then only 5 test vectors are used.

     If the matrix is not explicit, e.g., when estimating the norm of inv (A) given an LU factorization, onenormest applies A and its conjugate transpose through a pair of functions APPLY and APPLY_T, respectively, to a dense matrix of size N by T.  The implicit version requires an explicit dimension N.

     Returns the norm estimate EST, two vectors V and W related by norm (W, 1) = EST * norm (V, 1), and the number of iterations ITER.  The number of iterations is limited to 10 and is at least 2.

     References:

       * Nicholas J. Higham and Françoise Tisseur, "A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra." SIMAX vol 21, no 4, pp 1185-1201. http://dx.doi.org/10.1137/S0895479899356080

       * Nicholas J. Higham and Françoise Tisseur, "A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra." http://citeseer.ist.psu.edu/223007.html

     condest, norm, cond

 -- Function File: [EST, V] = condest (A, T)
 -- Function File: [EST, V] = condest (A, SOLVE, SOLVE_T, T)
 -- Function File: [EST, V] = condest (APPLY, APPLY_T, SOLVE, SOLVE_T, N, T)

     Estimate the 1-norm condition number of a matrix A using T test vectors using a randomized 1-norm estimator. If T exceeds 5, then only 5 test vectors are used.

     If the matrix is not explicit, e.g., when estimating the condition number of A given an LU factorization, condest uses the following functions:

     apply
          A*x for a matrix x of size N by T.

     apply_t
          A'*x for a matrix x of size N by T.

     solve
          A \ b for a matrix b of size N by T.

     solve_t
          A' \ b for a matrix b of size N by T.

     The implicit version requires an explicit dimension N.

     condest uses a randomized algorithm to approximate the 1-norms.

     condest returns the 1-norm condition estimate EST and a vector V satisfying norm (A*v, 1) == norm (A, 1) * norm (V, 1) / EST.  When EST is large, V is an approximate null vector.

     References:

       * Nicholas J. Higham and Françoise Tisseur, "A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra." SIMAX vol 21, no 4, pp 1185-1201. http://dx.doi.org/10.1137/S0895479899356080

       * Nicholas J. Higham and Françoise Tisseur, "A Block Algorithm for Matrix 1-Norm Estimation, with an Application to 1-Norm Pseudospectra." http://citeseer.ist.psu.edu/223007.html

     cond, norm, onenormest

 -- Loadable Function:   spparms ()
 -- Loadable Function: VALS = spparms ()
 -- Loadable Function: [KEYS, VALS] = spparms ()
 -- Loadable Function: VAL = spparms (KEY)
 -- Loadable Function:   spparms (VALS)
 -- Loadable Function:   spparms ('defaults')
 -- Loadable Function:   spparms ('tight')
 -- Loadable Function:   spparms (KEY, VAL)
     Sets or displays the parameters used by the sparse solvers and factorization functions.  The first four calls above get information about the current settings, while the others change the current settings.  The parameters are stored as pairs of keys and values, where the values are all floats and the keys are one of the following strings:

     spumoni
          Printing level of debugging information of the solvers (default 0)

     ths_rel
          Included for compatibility.  Not used.  (default 1)

     ths_abs
          Included for compatibility.  Not used.  (default 1)

     exact_d
          Included for compatibility.  Not used.  (default 0)

     supernd
          Included for compatibility.  Not used.  (default 3)

     rreduce
          Included for compatibility.  Not used.  (default 3)

     wh_frac
          Included for compatibility.  Not used.  (default 0.5)

     autommd
          Flag whether the LU/QR and the '\' and '/' operators will automatically use the sparsity preserving mmd functions (default 1)

     autoamd
          Flag whether the LU and the '\' and '/' operators will automatically use the sparsity preserving amd functions (default 1)

     piv_tol
          The pivot tolerance of the UMFPACK solvers (default 0.1)

     sym_tol
          The pivot tolerance of the UMFPACK symmetric solvers (default 0.001)

     bandden
          The density of non-zero elements in a banded matrix before it is treated by the LAPACK banded solvers (default 0.5)

     umfpack
          Flag whether the UMFPACK or mmd solvers are used for the LU, '\' and '/' operations (default 1)

     The value of individual keys can be set with spparms (KEY, VAL).  The default values can be restored with the special keyword 'defaults'.  The special keyword 'tight' can be used to set the mmd solvers to attempt for a sparser solution at the potential cost of longer running time.

 -- Loadable Function: P = sprank (S)

     Calculates the structural rank of a sparse matrix S.  Note that only the structure of the matrix is used in this calculation based on a Dulmage-Mendelsohn permutation to block triangular form.  As such the numerical rank of the matrix S is bounded by sprank (S) >= rank (S).  Ignoring floating point errors sprank (S) == rank (S). dmperm

 -- Loadable Function: [COUNT, H, PARENT, POST, R] = symbfact (S, TYP, MODE)

     Performs a symbolic factorization analysis on the sparse matrix S. Where

     S
          S is a complex or real sparse matrix.

     TYP
          Is the type of the factorization and can be one of

          sym
               Factorize S.  This is the default.

          col
               Factorize S' * S.

          row
               Factorize S * S'.

          lo
               Factorize S'

     MODE
          The default is to return the Cholesky factorization for R, and if MODE is 'L', the conjugate transpose of the Cholesky factorization is returned.  The conjugate transpose version is faster and uses less memory, but returns the same values for COUNT, H, PARENT and POST outputs.

     The output variables are

     COUNT
          The row counts of the Cholesky factorization as determined by TYP.

     H
          The height of the elimination tree.

     PARENT
          The elimination tree itself.

     POST
          A sparse boolean matrix whose structure is that of the Cholesky factorization as determined by TYP.

For non square matrices, the user can also utilize the spaugment function to find a least squares solution to a linear equation.

 -- Function File: S = spaugment (A, C)
     Creates the augmented matrix of A.  This is given by

          [C * eye(M, M),A; A', zeros(N,
          N)]

     This is related to the least squares solution of A \\ B, by

          S * [ R / C; x] = [B, zeros(N,
          columns(B)]

     where R is the residual error

          R = B - A * X

     As the matrix S is symmetric indefinite it can be factorized with lu, and the minimum norm solution can therefore be found without the need for a qr factorization.  As the residual error will be zeros (M, M) for under determined problems, and example can be

          m = 11; n = 10; mn = max(m ,n);
          a = spdiags ([ones(mn,1), 10*ones(mn,1), -ones(mn,1)],
                       [-1, 0, 1], m, n);
          x0 = a \ ones (m,1);
          s = spaugment (a);
          [L, U, P, Q] = lu (s);
          x1 = Q * (U \ (L \ (P  * [ones(m,1); zeros(n,1)])));
          x1 = x1(end - n + 1 : end);

     To find the solution of an overdetermined problem needs an estimate of the residual error R and so it is more complex to formulate a minimum norm solution using the spaugment function.

     In general the left division operator is more stable and faster than using the spaugment function.

Finally, the function eigs can be used to calculate a limited number of eigenvalues and eigenvectors based on a selection criteria and likewise for svds which calculates a limited number of singular values and vectors.

 -- Loadable Function: D = eigs (A)
 -- Loadable Function: D = eigs (A, K)
 -- Loadable Function: D = eigs (A, K, SIGMA)
 -- Loadable Function: D = eigs (A, K, SIGMA,OPTS)
 -- Loadable Function: D = eigs (A, B)
 -- Loadable Function: D = eigs (A, B, K)
 -- Loadable Function: D = eigs (A, B, K, SIGMA)
 -- Loadable Function: D = eigs (A, B, K, SIGMA, OPTS)
 -- Loadable Function: D = eigs (AF, N)
 -- Loadable Function: D = eigs (AF, N, B)
 -- Loadable Function: D = eigs (AF, N, K)
 -- Loadable Function: D = eigs (AF, N, B, K)
 -- Loadable Function: D = eigs (AF, N, K, SIGMA)
 -- Loadable Function: D = eigs (AF, N, B, K, SIGMA)
 -- Loadable Function: D = eigs (AF, N, K, SIGMA, OPTS)
 -- Loadable Function: D = eigs (AF, N, B, K, SIGMA, OPTS)
 -- Loadable Function: [V, D] = eigs (A, ...)
 -- Loadable Function: [V, D] = eigs (AF, N, ...)
 -- Loadable Function: [V, D, FLAG] = eigs (A, ...)
 -- Loadable Function: [V, D, FLAG] = eigs (AF, N, ...)
     Calculate a limited number of eigenvalues and eigenvectors of A, based on a selection criteria.  The number eigenvalues and eigenvectors to calculate is given by K whose default value is 6.

     By default eigs solve the equation A * v = lambda * v , where lambda is a scalar representing one of the eigenvalues, and v is the corresponding eigenvector.  If given the positive definite matrix B then eigs solves the general eigenvalue equation A * v = lambda * B * v .

     The argument SIGMA determines which eigenvalues are returned. SIGMA can be either a scalar or a string.  When SIGMA is a scalar, the K eigenvalues closest to SIGMA are returned.  If SIGMA is a string, it must have one of the values

     'lm'
          Largest magnitude (default).

     'sm'
          Smallest magnitude.

     'la'
          Largest Algebraic (valid only for real symmetric problems).

     'sa'
          Smallest Algebraic (valid only for real symmetric problems).

     'be'
          Both ends, with one more from the high-end if K is odd (valid only for real symmetric problems).

     'lr'
          Largest real part (valid only for complex or unsymmetric problems).

     'sr'
          Smallest real part (valid only for complex or unsymmetric problems).

     'li'
          Largest imaginary part (valid only for complex or unsymmetric problems).

     'si'
          Smallest imaginary part (valid only for complex or unsymmetric problems).

     If OPTS is given, it is a structure defining some of the options that eigs should use.  The fields of the structure OPTS are

     issym
          If AF is given, then flags whether the function AF defines a symmetric problem.  It is ignored if A is given.  The default is false.

     isreal
          If AF is given, then flags whether the function AF defines a real problem.  It is ignored if A is given.  The default is true.

     tol
          Defines the required convergence tolerance, given as tol * norm (A). The default is eps.

     maxit
          The maximum number of iterations.  The default is 300.

     p
          The number of Lanzcos basis vectors to use.  More vectors will result in faster convergence, but a larger amount of memory.  The optimal value of 'p' is problem dependent and should be in the range K to N.  The default value is 2 * K.

     v0
          The starting vector for the computation.  The default is to have ARPACK randomly generate a starting vector.

     disp
          The level of diagnostic printout.  If disp is 0 then there is no printout.  The default value is 1.

     cholB
          Flag if chol (B) is passed rather than B.  The default is false.

     permB
          The permutation vector of the Cholesky factorization of B if cholB is true.  That is chol ( B (permB, permB)).  The default is 1:N.

     It is also possible to represent A by a function denoted AF. AF must be followed by a scalar argument N defining the length of the vector argument accepted by AF.  AF can be passed either as an inline function, function handle or as a string.  In the case where AF is passed as a string, the name of the string defines the function to use.

     AF is a function of the form function y = af (x), y = ...; endfunction, where the required return value of AF is determined by the value of SIGMA, and are

     A * x
          If SIGMA is not given or is a string other than 'sm'.

     A \ x
          If SIGMA is 'sm'.

     (A - sigma * I) \ x
          for standard eigenvalue problem, where I is the identity matrix of the same size as A.  If SIGMA is zero, this reduces the A \ x.

     (A - sigma * B) \ x
          for the general eigenvalue problem.

     The return arguments of eigs depends on the number of return arguments.  With a single return argument, a vector D of length K is returned, represent the K eigenvalues that have been found.  With two return arguments, V is a N-by-K matrix whose columns are the K eigenvectors corresponding to the returned eigenvalues.  The eigenvalues themselves are then returned in D in the form of a N-by-K matrix, where the elements on the diagonal are the eigenvalues.

     Given a third return argument FLAG, eigs also returns the status of the convergence.  If FLAG is 0, then all eigenvalues have converged, otherwise not.

     This function is based on the ARPACK package, written by R Lehoucq, K Maschhoff, D Sorensen and C Yang.  For more information see http://www.caam.rice.edu/software/ARPACK/.

eig, svds

 -- Function File: S = svds (A)
 -- Function File: S = svds (A, K)
 -- Function File: S = svds (A, K, SIGMA)
 -- Function File: S = svds (A, K, SIGMA, OPTS)
 -- Function File: [U, S, V, FLAG] = svds (...)

     Find a few singular values of the matrix A.  The singular values are calculated using

          [M, N] = size(A)
          S = eigs([sparse(M, M), A; ...
                          A', sparse(N, N)])

     The eigenvalues returned by eigs correspond to the singular values of A.  The number of singular values to calculate is given by K, whose default value is 6.

     The argument SIGMA can be used to specify which singular values to find.  SIGMA can be either the string 'L', the default, in which case the largest singular values of A are found.  Otherwise SIGMA should be a real scalar, in which case the singular values closest to SIGMA are found.  Note that for relatively small values of SIGMA, there is the chance that the requested number of singular values are not returned.  In that case SIGMA should be increased.

     If OPTS is given, then it is a structure that defines options that svds will pass to EIGS.  The possible fields of this structure are therefore determined by eigs.  By default three fields of this structure are set by svds.

     tol
          The required convergence tolerance for the singular values.  eigs is passed TOL divided by sqrt(2).  The default value is 1e-10.

     maxit
          The maximum number of iterations.  The default is 300.

     disp
          The level of diagnostic printout.  If disp is 0 then there is no printout.  The default value is 0.

     If more than one output argument is given, then svds also calculates the left and right singular vectors of A.  FLAG is used to signal the convergence of svds.  If svds converges to the desired tolerance, then FLAG given by

          norm (A * V - U * S, 1) <= ...
                  TOL * norm (A, 1)

     will be zero.

eigs

Iterative Techniques applied to sparse matrices
===============================================

The left division \ and right division / operators, discussed in the previous section, use direct solvers to resolve a linear equation of the form X = A \ B or X = B / A.  Octave equally includes a number of functions to solve sparse linear equations using iterative techniques.

 -- Function File: X = pcg (A, B, TOL, MAXIT, M1, M2, X0, ...)
 -- Function File: [X, FLAG, RELRES, ITER, RESVEC, EIGEST] = pcg (...)

     Solves the linear system of equations A * X = B by means of the Preconditioned Conjugate Gradient iterative method.  The input arguments are

       * A can be either a square (preferably sparse) matrix or a function handle, inline function or string containing the name of a function which computes A * X.  In principle A should be symmetric and positive definite; if pcg finds A to not be positive definite, you will get a warning message and the FLAG output parameter will be set.

       * B is the right hand side vector.

       * TOL is the required relative tolerance for the residual error, B - A * X.  The iteration stops if norm (B - A * X) <= TOL * norm (B - A * X0).  If TOL is empty or is omitted, the function sets TOL = 1e-6 by default.

       * MAXIT is the maximum allowable number of iterations; if [] is supplied for maxit, or pcg has less arguments, a default value equal to 20 is used.

       * M = M1 * M2 is the (left) preconditioning matrix, so that the iteration is (theoretically) equivalent to solving by pcg P * X = M \ B, with P = M \ A. Note that a proper choice of the preconditioner may dramatically improve the overall performance of the method.  Instead of matrices M1 and M2, the user may pass two functions which return the results of applying the inverse of M1 and M2 to a vector (usually this is the preferred way of using the preconditioner). If [] is supplied for M1, or M1 is omitted, no preconditioning is applied.  If M2 is omitted, M = M1 will be used as preconditioner.

       * X0 is the initial guess.  If X0 is empty or omitted, the function sets X0 to a zero vector by default.

     The arguments which follow X0 are treated as parameters, and passed in a proper way to any of the functions (A or M) which are passed to pcg.  See the examples below for further details.  The output arguments are

       * X is the computed approximation to the solution of A * X = B.

       * FLAG reports on the convergence.  FLAG = 0 means the solution converged and the tolerance criterion given by TOL is satisfied.  FLAG = 1 means that the MAXIT limit for the iteration count was reached.  FLAG = 3 reports that the (preconditioned) matrix was found not positive definite.

       * RELRES is the ratio of the final residual to its initial value, measured in the Euclidean norm.

       * ITER is the actual number of iterations performed.

       * RESVEC describes the convergence history of the method. RESVEC (i,1) is the Euclidean norm of the residual, and RESVEC (i,2) is the preconditioned residual norm, after the (I-1)-th iteration, I = 1, 2, ..., ITER+1.  The preconditioned residual norm is defined as norm (R) ^ 2 = R' * (M \ R) where R = B - A * X, see also the description of M.  If EIGEST is not required, only RESVEC (:,1) is returned.

       * EIGEST returns the estimate for the smallest EIGEST (1) and largest EIGEST (2) eigenvalues of the preconditioned matrix P = M \ A.  In particular, if no preconditioning is used, the estimates for the extreme eigenvalues of A are returned.  EIGEST (1) is an overestimate and EIGEST (2) is an underestimate, so that EIGEST (2) / EIGEST (1) is a lower bound for cond (P, 2), which nevertheless in the limit should theoretically be equal to the actual value of the condition number. The method which computes EIGEST works only for symmetric positive definite A and M, and the user is responsible for verifying this assumption.

     Let us consider a trivial problem with a diagonal matrix (we exploit the sparsity of A)

          	n = 10;
          	a = diag (sparse (1:n));
          	b = rand (n, 1);
               [l, u, p, q] = luinc (a, 1.e-3);

     EXAMPLE 1: Simplest use of pcg

            x = pcg(A,b)

     EXAMPLE 2: pcg with a function which computes A * X

            function y = apply_a (x)
              y = [1:N]'.*x;
            endfunction

            x = pcg ("apply_a", b)

     EXAMPLE 3: pcg with a preconditioner: L * U

          x = pcg (a, b, 1.e-6, 500, l*u);

     EXAMPLE 4: pcg with a preconditioner: L * U. Faster than EXAMPLE 3 since lower and upper triangular matrices are easier to invert

          x = pcg (a, b, 1.e-6, 500, l, u);

     EXAMPLE 5: Preconditioned iteration, with full diagnostics.  The preconditioner (quite strange, because even the original matrix A is trivial) is defined as a function

            function y = apply_m (x)
              k = floor (length (x) - 2);
              y = x;
              y(1:k) = x(1:k)./[1:k]';
            endfunction

            [x, flag, relres, iter, resvec, eigest] = ...
                               pcg (a, b, [], [], "apply_m");
            semilogy (1:iter+1, resvec);

     EXAMPLE 6: Finally, a preconditioner which depends on a parameter K.

            function y = apply_M (x, varargin)
            K = varargin{1};
            y = x;
            y(1:K) = x(1:K)./[1:K]';
            endfunction

            [x, flag, relres, iter, resvec, eigest] = ...
                 pcg (A, b, [], [], "apply_m", [], [], 3)

     REFERENCES

     [1] C.T.Kelley, 'Iterative methods for linear and nonlinear equations', SIAM, 1995 (the base PCG algorithm)

     [2] Y.Saad, 'Iterative methods for sparse linear systems', PWS 1996 (condition number estimate from PCG) Revised version of this book is available online at http://www-users.cs.umn.edu/~saad/books.html

     sparse, pcr

 -- Function File: X = pcr (A, B, TOL, MAXIT, M, X0, ...)
 -- Function File: [X, FLAG, RELRES, ITER, RESVEC] = pcr (...)

     Solves the linear system of equations A * X = B by means of the Preconditioned Conjugate Residuals iterative method.  The input arguments are

       * A can be either a square (preferably sparse) matrix or a function handle, inline function or string containing the name of a function which computes A * X.  In principle A should be symmetric and non-singular; if pcr finds A to be numerically singular, you will get a warning message and the FLAG output parameter will be set.

       * B is the right hand side vector.

       * TOL is the required relative tolerance for the residual error, B - A * X.  The iteration stops if norm (B - A * X) <= TOL * norm (B - A * X0).  If TOL is empty or is omitted, the function sets TOL = 1e-6 by default.

       * MAXIT is the maximum allowable number of iterations; if [] is supplied for maxit, or pcr has less arguments, a default value equal to 20 is used.

       * M is the (left) preconditioning matrix, so that the iteration is (theoretically) equivalent to solving by pcr P * X = M \ B, with P = M \ A. Note that a proper choice of the preconditioner may dramatically improve the overall performance of the method.  Instead of matrix M, the user may pass a function which returns the results of applying the inverse of M to a vector (usually this is the preferred way of using the preconditioner).  If [] is supplied for M, or M is omitted, no preconditioning is applied.

       * X0 is the initial guess.  If X0 is empty or omitted, the function sets X0 to a zero vector by default.

     The arguments which follow X0 are treated as parameters, and passed in a proper way to any of the functions (A or M) which are passed to pcr.  See the examples below for further details.  The output arguments are

       * X is the computed approximation to the solution of A * X = B.

       * FLAG reports on the convergence.  FLAG = 0 means the solution converged and the tolerance criterion given by TOL is satisfied.  FLAG = 1 means that the MAXIT limit for the iteration count was reached.  FLAG = 3 reports t pcr breakdown, see [1] for details.

       * RELRES is the ratio of the final residual to its initial value, measured in the Euclidean norm.

       * ITER is the actual number of iterations performed.

       * RESVEC describes the convergence history of the method, so that RESVEC (i) contains the Euclidean norms of the residual after the (I-1)-th iteration, I = 1,2, ..., ITER+1.

     Let us consider a trivial problem with a diagonal matrix (we exploit the sparsity of A)

          	n = 10;
          	a = sparse (diag (1:n));
          	b = rand (N, 1);

     EXAMPLE 1: Simplest use of pcr

            x = pcr(A, b)

     EXAMPLE 2: pcr with a function which computes A * X.

            function y = apply_a (x)
              y = [1:10]'.*x;
            endfunction

            x = pcr ("apply_a", b)

     EXAMPLE 3:  Preconditioned iteration, with full diagnostics.  The preconditioner (quite strange, because even the original matrix A is trivial) is defined as a function

            function y = apply_m (x)
              k = floor (length(x)-2);
              y = x;
              y(1:k) = x(1:k)./[1:k]';
            endfunction

            [x, flag, relres, iter, resvec] = ...
                               pcr (a, b, [], [], "apply_m")
            semilogy([1:iter+1], resvec);

     EXAMPLE 4: Finally, a preconditioner which depends on a parameter K.

            function y = apply_m (x, varargin)
              k = varargin{1};
              y = x; y(1:k) = x(1:k)./[1:k]';
            endfunction

            [x, flag, relres, iter, resvec] = ...
                               pcr (a, b, [], [], "apply_m"', [], 3)

     REFERENCES

     [1] W. Hackbusch, "Iterative Solution of Large Sparse Systems of Equations", section 9.5.4; Springer, 1994

     sparse, pcg

The speed with which an iterative solver converges to a solution can be accelerated with the use of a pre-conditioning matrix M.  In this case the linear equation M^-1 * X = M^-1 * A \ B is solved instead.  Typical pre-conditioning matrices are partial factorizations of the original matrix.

 -- Loadable Function: [L, U, P, Q] = luinc (A, '0')
 -- Loadable Function: [L, U, P, Q] = luinc (A, DROPTOL)
 -- Loadable Function: [L, U, P, Q] = luinc (A, OPTS)
     Produce the incomplete LU factorization of the sparse matrix A. Two types of incomplete factorization are possible, and the type is determined by the second argument to luinc.

     Called with a second argument of '0', the zero-level incomplete LU factorization is produced.  This creates a factorization of A where the position of the non-zero arguments correspond to the same positions as in the matrix A.

     Alternatively, the fill-in of the incomplete LU factorization can be controlled through the variable DROPTOL or the structure OPTS.  The UMFPACK multifrontal factorization code by Tim A. Davis is used for the incomplete LU factorization, (availability http://www.cise.ufl.edu/research/sparse/umfpack/)

     DROPTOL determines the values below which the values in the LU factorization are dropped and replaced by zero.  It must be a positive scalar, and any values in the factorization whose absolute value are less than this value are dropped, expect if leaving them increase the sparsity of the matrix.  Setting DROPTOL to zero results in a complete LU factorization which is the default.

     OPTS is a structure containing one or more of the fields

     droptol
          The drop tolerance as above.  If OPTS only contains droptol then this is equivalent to using the variable DROPTOL.

     milu
          A logical variable flagging whether to use the modified incomplete LU factorization.  In the case that milu is true, the dropped values are subtracted from the diagonal of the matrix U of the factorization. The default is false.

     udiag
          A logical variable that flags whether zero elements on the diagonal of U should be replaced with DROPTOL to attempt to avoid singular factors.  The default is false.

     thresh
          Defines the pivot threshold in the interval [0,1].  Values outside that range are ignored.

     All other fields in OPTS are ignored.  The outputs from luinc are the same as for lu.

     Given the string argument 'vector', luinc returns the values of P Q as vector values. sparse, lu

Real Life Example of the use of Sparse Matrices
===============================================

A common application for sparse matrices is in the solution of Finite Element Models.  Finite element models allow numerical solution of partial differential equations that do not have closed form solutions, typically because of the complex shape of the domain.

In order to motivate this application, we consider the boundary value Laplace equation.  This system can model scalar potential fields, such as heat or electrical potential.  Given a medium Omega with boundary dOmega . At all points on the dOmega the boundary conditions are known, and we wish to calculate the potential in Omega . Boundary conditions may specify the potential (Dirichlet boundary condition), its normal derivative across the boundary (Neumann boundary condition), or a weighted sum of the potential and its derivative (Cauchy boundary condition).

In a thermal model, we want to calculate the temperature in Omega and know the boundary temperature (Dirichlet condition) or heat flux (from which we can calculate the Neumann condition by dividing by the thermal conductivity at the boundary).  Similarly, in an electrical model, we want to calculate the voltage in Omega and know the boundary voltage (Dirichlet) or current (Neumann condition after diving by the electrical conductivity). In an electrical model, it is common for much of the boundary to be electrically isolated; this is a Neumann boundary condition with the current equal to zero.

The simplest finite element models will divide Omega into simplexes (triangles in 2D, pyramids in 3D).

The following example creates a simple rectangular 2D electrically conductive medium with 10 V and 20 V imposed on opposite sides (Dirichlet boundary conditions).  All other edges are electrically isolated.

        node_y= [1;1.2;1.5;1.8;2]*ones(1,11);
        node_x= ones(5,1)*[1,1.05,1.1,1.2, ...
                  1.3,1.5,1.7,1.8,1.9,1.95,2];
        nodes= [node_x(:), node_y(:)];

        [h,w]= size(node_x);
        elems= [];
        for idx= 1:w-1
          widx= (idx-1)*h;
          elems= [elems; ...
            widx+[(1:h-1);(2:h);h+(1:h-1)]'; ...
            widx+[(2:h);h+(2:h);h+(1:h-1)]' ];
        endfor

        E= size(elems,1); # No. of simplices
        N= size(nodes,1); # No. of vertices
        D= size(elems,2); # dimensions+1

This creates a N-by-2 matrix nodes and a E-by-3 matrix elems with values, which define finite element triangles:

       nodes(1:7,:)'
         1.00 1.00 1.00 1.00 1.00 1.05 1.05 ...
         1.00 1.20 1.50 1.80 2.00 1.00 1.20 ...

       elems(1:7,:)'
         1    2    3    4    2    3    4 ...
         2    3    4    5    7    8    9 ...
         6    7    8    9    6    7    8 ...

Using a first order FEM, we approximate the electrical conductivity distribution in Omega as constant on each simplex (represented by the vector conductivity). Based on the finite element geometry, we first calculate a system (or stiffness) matrix for each simplex (represented as 3-by-3 elements on the diagonal of the element-wise system matrix SE.  Based on SE and a N-by-DE connectivity matrix C, representing the connections between simplices and vertices, the global connectivity matrix S is calculated.

       # Element conductivity
       conductivity= [1*ones(1,16), ...
              2*ones(1,48), 1*ones(1,16)];

       # Connectivity matrix
       C = sparse ((1:D*E), reshape (elems', ...
              D*E, 1), 1, D*E, N);

       # Calculate system matrix
       Siidx = floor ([0:D*E-1]'/D) * D * ...
              ones(1,D) + ones(D*E,1)*(1:D) ;
       Sjidx = [1:D*E]'*ones(1,D);
       Sdata = zeros(D*E,D);
       dfact = factorial(D-1);
       for j=1:E
          a = inv([ones(D,1), ...
              nodes(elems(j,:), :)]);
          const = conductivity(j) * 2 / ...
              dfact / abs(det(a));
          Sdata(D*(j-1)+(1:D),:) = const * ...
              a(2:D,:)' * a(2:D,:);
       endfor
       # Element-wise system matrix
       SE= sparse(Siidx,Sjidx,Sdata);
       # Global system matrix
       S= C'* SE *C;

The system matrix acts like the conductivity S in Ohm's law S * V = I. Based on the Dirichlet and Neumann boundary conditions, we are able to solve for the voltages at each vertex V.

       # Dirichlet boundary conditions
       D_nodes=[1:5, 51:55];
       D_value=[10*ones(1,5), 20*ones(1,5)];

       V= zeros(N,1);
       V(D_nodes) = D_value;
       idx = 1:N; # vertices without Dirichlet
                  # boundary condns
       idx(D_nodes) = [];

       # Neumann boundary conditions.  Note that
       # N_value must be normalized by the
       # boundary length and element conductivity
       N_nodes=[];
       N_value=[];

       Q = zeros(N,1);
       Q(N_nodes) = N_value;

       V(idx) = S(idx,idx) \ ( Q(idx) - ...
                 S(idx,D_nodes) * V(D_nodes));

Finally, in order to display the solution, we show each solved voltage value in the z-axis for each simplex vertex.

       elemx = elems(:,[1,2,3,1])';
       xelems = reshape (nodes(elemx, 1), 4, E);
       yelems = reshape (nodes(elemx, 2), 4, E);
       velems = reshape (V(elemx), 4, E);
       plot3 (xelems,yelems,velems,'k');
       print ('grid.eps');