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

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  (Y. 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.  It is similar to the Yale format. (http://en.wikipedia.org/wiki/Sparse_matrix#Yale_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.

(spdiags)

(speye)

(spones)

(sprand)

(sprandn)

(sprandsym)

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,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 pre-assign 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.

(full)

(spalloc)

(sparse)

(spconvert)

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.  See External Code Interface, for a full description of the techniques involved.

Finding 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.

(issparse)

(nnz)

(nonzeros)

(nzmax)

(spstats)

When solving linear equations involving sparse matrices Octave determines the means to solve the equation based on the type of the matrix (see 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

shows 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.

(spy)

(etree)

(etreeplot)

(gplot)

(treeplot)

(treelayout)

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

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

Many Octave functions have been overloaded to work with either sparse or full matrices.  There is no difference in calling convention when using an overloaded function with a sparse matrix, however, there is also no access to potentially sparse-specific features.  At any time the sparse matrix specific version of a 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 in the rest of this article, the specific sparse versions of 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.

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.

(sparse_auto_mutate)

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.

First, 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);. See fig:simplechol. 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, p, q] = chol (A); spy (r).

Figure: Structure of the unpermuted Cholesky factorization of the above matrix.

Figure: Structure of the permuted Cholesky factorization of the above matrix.

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.

(amd)

(ccolamd)

(colamd)

(colperm)

(csymamd)

(dmperm)

(symamd)

(symrcm)

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 an 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 band divided by the total number of values in the full band.  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 decomposition, 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.

(normest)

(onenormest)

(condest)

(spparms)

(sprank)

(symbfact)

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

(spaugment)

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.

(eigs)

(svds)

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.

(pcg)

(pcr)

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.

(luinc)

Real Life Example using 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 2-D 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";