@c Copyright (C) 2004, 2005 David Bateman
@c This is part of the Octave manual.
@c For copying conditions, see the file gpl.texi.

@node Sparse Matrices 
@chapter Sparse Matrices

@menu
* Basics:: The Creation and Manipulation of Sparse Matrices
* Graph Theory:: Graphs are their use with Sparse Matrices
* Sparse Linear Algebra:: Linear Algebra on Sparse Matrices
* Iterative Techniques:: Iterative Techniques applied to Sparse Matrices
* Oct-Files:: Using Sparse Matrices in Oct-files
* License:: Licensing of Third Party Software
* Function Reference:: Documentation from the Specific Sparse Functions
@end menu

@node Basics, Graph Theory, Sparse Matrices, Sparse Matrices
@section 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 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 done 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.

THIS DOCUMENT STILL HAS LARGE BLANKS. PLEASE FILL THEM IN. LOOK FOR
THE TAG "WRITE ME"

@menu
* Storage:: Storage of Sparse Matrices
* Creation:: Creating Sparse Matrices
* Operators and Functions:: Basic Operators and Functions on Sparse Matrices
* Information:: Finding out Information about Sparse Matrices
@end menu

@node Storage, Creation, Basics, Basics
@subsection 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 compelxity
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 @footnote{Youcef Saad "SPARSKIT: A basic toolkit
for sparse matrix computation", 1994,
@url{ftp://ftp.cs.umn.edu/dept/sparse/SPARSKIT2/DOC/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 the of the 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 @sc{Octave} is 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 simplication in the code, in that their is no special case
for the first or last columns. A short example, demonstrating this in
C is.

@example
  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));
@end example

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

@example
@group
    1   2   0  0
    0   0   0  3
    0   0   0  4
@end group
@end example

The non-zero elements of this matrix are

@example
@group
   (1, 1)  @result{} 1
   (1, 2)  @result{} 2
   (2, 4)  @result{} 3
   (3, 4)  @result{} 4
@end group
@end example

This will be stored as three vectors @var{cidx}, @var{ridx} and @var{data},
representing the column indexing, row indexing and data respectively. The
contents of these three vectors for the above matrix will be

@example
@group
  @var{cidx} = [0, 2, 2, 4]
  @var{ridx} = [0, 0, 1, 2]
  @var{data} = [1, 2, 3, 4]
@end group
@end example

Note that this is the representation of these elements with the first row
and column assumed to start as zero. Thus the number of elements in the 
@var{i}-th column is given by @code{@var{cidx} (@var{i} + 1) - 
@var{cidx} (@var{i})}.

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. @sc{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 @sc{Octave} is that 
all elements in the rows are stored in increasing order of their row
index. This makes certain operations, later be faster. However, it imposes
the need to sort the elements on the creation of sparse matrices. Having
dis-ordered 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.

@node Creation, Operators and Functions, Storage, Basics
@subsection Creating Sparse Matrices

There are several means to create sparse matrix.

@table @asis
@item Returned from a function
There are many functions that directly return sparse matrices. These include
@dfn{speye}, @dfn{sprand}, @dfn{spdiag}, etc.
@item Constructed from matrices or vectors
The function @dfn{sparse} allows a sparse matrix to be constructed from 
three vectors representing the row, column and data. ALternatively, the
function @dfn{spconvert} uses a three column matrix format to allow easy
importation of data from elsewhere.
@item Created and then filled
The function @dfn{sparse} or @dfn{spalloc} can be used to create an empty
matrix that is then filled by the user
@item From a user binary program
The user can directly create the sparse matrix within an oct-file.
@end table

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
@code{speye (@var{n})} or @code{speye (@var{r}, @var{c})}, which
creates an @var{n}-by-@var{n} or @var{r}-by-@var{c} sparse identity
matrix.

Another typical sparse matrix that is often needed is a random distribution
of random elements. The functions @dfn{sprand} and @dfn{sprandn} perform
this for uniform and random distributions of elements. They have exactly
the same calling convention as @code{sprand (@var{r}, @var{c}, @var{d})},
which creates an @var{r}-by-@var{c} sparse matrix a density of filled
elements of @var{d}.

Other functions of interest that directly creates a sparse matrix, are
@dfn{spdiag} or its generalization @dfn{spdiags}, that can take the
definition of the diagonals of the matrix and create the sparse matrix 
that corresponds to this. For example

@c FIXME, when spdiag is overloaded to diag the below won't work. 

@example
s = spdiag (randn(1,n), -1);
@end example

creates a sparse (@var{n}+1)-by-(@var{n}+1) sparse matrix with a single
diagonal defined.

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

@example
 function x = foo (r, j)
  idx = randperm (r);
  x = ([zeros(r-2,1); rand(2,1)])(idx);
 endfunction

 ri = [];
 ci = [];
 d = [];

 for j=1:c
    dtmp = foo (r, j);  # Returns vector of length r with (:,j) values
    idx = find (dtmp != 0.);
    ri = [ri; idx];
    ci = [ci; j*ones(length(idx),1)]; 
    d = [d; dtmp(idx)];
 endfor
 s = sparse (ri, ci, d, r, c);
@end example

creates an @var{r}-by-@var{c} sparse matrix with a random distribution of
2 elements per row. The elements of the vectors do not need to be sorted in
any particular order as @sc{Octave} will store them prior to storing the
data. However, per sorting the data will make teh creation of the sparse
matrix faster.

The function @dfn{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

@example
s = spconvert ([1 2 3 4; 1 3 4 4; 1 2 3 0]')
@result{} Compressed Column Sparse (rows=4, cols=4, nnz=3)
      (1 , 1) -> 1
      (2 , 3) -> 2
      (3 , 4) -> 3
@end example

An example of creating and filling a matrix might be

@example
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
@end example

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

The above problem can be avoided in oct-files. However, the
construction of a sparse matrix from an oct-file is more complex than
can be discussed in this brief introduction, and you are referred to
section @ref{Oct-Files}, to have a full description of the techniques
involved.

@node Operators and Functions, Information, Creation, Basics
@subsection Basic Operators and Functions on Sparse Matrices

@menu
* Functions:: Operators and Functions
* ReturnType:: The Return Types of Operators and Functions
* MathConsiderations:: Mathematical Considerations
@end menu

@node Functions, ReturnType, Operators and Functions, Operators and Functions
@subsubsection Operators and Functions

WRITE ME

@node ReturnType, MathConsiderations, Functions, Operators and Functions
@subsubsection The Return Types of Operators and Functions

The two basic reason 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 numberof 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

@example
speye(3) + 0
@result{}   1  0  0
  0  1  0
  0  0  1
@end example

returns a full matrix as can be seen. Additionally 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 @code{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 as rare, but can be artificially created, for
example @code{(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 behaviour of @sc{Octave}'s
sparse matrices will cause a problem. That is in the handling of the
@dfn{diag} function. Whether @dfn{diag} returns a sparse or full matrix
depending on the type of its input arguments. So 

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

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

@node MathConsiderations, , ReturnType, Operators and Functions
@subsubsection 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

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

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

Likewise @code{@var{s} .^ -2} involves terms terms like @code{0 .^ -2} which
is infinity, and so @code{@var{s} .^ -2} is equally a full matrix.

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

The above behaviour 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

@example
 a = 0 ./ [-1, 1; 1, -1];
 b = 1 ./ a
 @result{} -Inf            Inf
     Inf           -Inf
 c = 1 ./ sparse (a)
 @result{}  Inf            Inf
     Inf            Inf
@end example
 
To correct this behaviour 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
efficient, and so the user is warned that calculations where the sign-bit
of zero is important must not be done using sparse matrices.

Also discuss issues of fill-in. Discuss symamd etc, and mention randperm
that is included  elsewhere in the docs...

WRITE ME

@node Information, , Operators and Functions, Basics
@subsection Finding out Information about Sparse Matrices

Talk about the spy, spstats, nnz, spparms, etc function

WRITE ME

@node Graph Theory, Sparse Linear Algebra, Basics, Sparse Matrices
@section Graphs are their use with Sparse Matrices

Someone who knows more about this than me should write this...

WRITE ME

@node Sparse Linear Algebra, Iterative Techniques, Graph Theory, Sparse Matrices
@section Linear Algebra on Sparse Matrices

@sc{Octave} includes a poly-morphic solver for sparse matrices, where 
the exact solver used to factorize the matrix, depends on the properties
of the sparse matrix, itself. As the cost of determining the matrix type
is small relative to the cost of factorizing the matrix itself, the matrix
type is re-determined each time it is used in a linear equation.

The selection tree for how the linear equation is solve is

@enumerate 1
@item If the matrix is not square go to 9.

@item If the matrix is diagonal, solve directly and goto 9

@item If the matrix is a permuted diagonal, solve directly taking into
account the permutations. Goto 9

@item If the matrix is banded and if the band density is less than that
given by @code{spparms ("bandden")} continue, else goto 5.

@enumerate a
@item If the matrix is tridiagonal and the right-hand side is not sparse 
continue, else goto 4b.

@enumerate
@item If the matrix is hermitian, with a positive real diagonal, attempt
      Cholesky factorization using @sc{Lapack} xPTSV.

@item If the above failed or the matrix is not hermitian with a positive
      real diagonal use Gaussian elimination with pivoting using 
      @sc{Lapack} xGTSV, and goto 9.
@end enumerate

@item If the matrix is hermitian with a positive real diagonal, attempt
      Cholesky factorization using @sc{Lapack} xPBTRF.

@item if the above failed or the matrix is not hermitian with a positive
      real diagonal use Gaussian elimination with pivoting using 
      @sc{Lapack} xGBTRF, and goto 9.
@end enumerate

@item If the matrix is upper or lower triangular perform a sparse forward
or backward subsitution, and goto 9

@item If the matrix is a permuted upper or lower triangular matrix, perform
a sparse forward or backward subsitution, and goto 9

FIXME: Detection of permuted triangular matrices not written yet, and so
       the code itself is not tested either

@item If the matrix is hermitian with a real positive diagonal, attempt
sparse Cholesky factorization.

FIXME: Detection of positive definite matrices written and tested, but 
   Cholesky factorization isn't yet written

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

@item If the matrix is not square, or any of the previous solvers flags
a singular or near singular matrix, find a minimum norm solution

FIXME: QR solvers not yet written.

@end enumerate

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 @dfn{spparms} to set @code{bandden}
to 1 (i.e. @code{spparms ("bandden", 1)}).

All of the solvers above, expect the banded solvers, 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 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 for banded matrices, and therefore unless the factorization
is exactly singular, these numerical instabilities won't be detected.
In cases where, this might be a problem the user is recommended to disable
the banded solvers as above, at a significant cost in terms of speed.

@node Iterative Techniques, Oct-Files, Sparse Linear Algebra, Sparse Matrices
@section Iterative Techniques applied to sparse matrices

WRITE ME

@node Oct-Files, License, Iterative Techniques, Sparse Matrices
@section Using Sparse Matrices in Oct-files

An oct-file is a means of writing an @sc{Octave} function in a compilable
language like C++, rather than as a script file. This results in a
significant acceleration in the code.  It is not the purpose of this
section to discuss how to write an oct-file, or discuss what they
are. There are already two @footnote{Paul Thomas "Dal Segno al Coda 
- The octave dynamically linked function cookbook", 
@url{http://perso.wanadoo.fr/prthomas/intro.html}, and Cristophe Spiel 
"Del Coda Al Fine - Pushing Octave's Limits", 
@url{http://octave.sourceforge.net/coda/coda.pdf}} very good
references on oct-files themselves. Users who are not familiar with
oct-files are urged to read these references to fully understand this
chapter. The examples discussed here assume that the oct-file is written 
entirely in C++.

There are three classes of sparse objects that are of interest to the
user.

@table @asis
@item SparseMatrix
A double precision sparse matrix class
@item SparseComplexMatrix
A Complex sparse matrix class
@item SparseBoolMatrix
A boolen sparse matrix class
@end table

All of these classes inherit from the @code{Sparse<T>} template class,
and so all have similar capabilities and usage. The @code{Sparse<T>}
class was based on @sc{Octave} @code{Array<T>} class, and so users familar
with @sc{Octave}'s Array classes will be comfortable with the use of
the sparse classes.

The sparse classes will not be entirely described in this section, due
to their similar with the existing Array classes. However, there are a
few differences due the different nature of sparse objects, and these
will be described. Firstly, although it is fundamentally possible to
have N-dimensional sparse objects, the @sc{Octave} sparse classes do
not allow them at this time. So all operations of the sparse classes
must be 2-dimensional.  This means that in fact @code{SparseMatrix} is
similar to @sc{Octave}'s @code{Matrix} class rather than its
@code{NDArray} class.

@menu
* OctDifferences:: The Differences between the Array and Sparse Classes
* OctCreation:: Creating Spare Matrices in Oct-Files
* OctUse:: Using Sparse Matrices in Oct-Files
@end menu

@node OctDifferences, OctCreation, Oct-Files, Oct-Files
@subsection The Differences between the Array and Sparse Classes

The number of elements in a sparse matrix is considered to be the number
of non-zero elements rather than the product of the dimensions. Therefore

@example
  SparseMatrix sm;
  @dots{}
  int nel = sm.nelem ();
@end example

returns the number of non-zero elements. If the user really requires the 
number of elements in the matrix, including the non-zero elements, they
should use @code{numel} rather than @code{nelem}. Note that for very 
large matrices, where the product of the two dimensions is large that
the representation of the an unsigned int, then @code{numel} can overflow.
An example is @code{speye(1e6)} which will create a matrix with a million
rows and columns, but only a million non-zero elements. Therefore the
number of rows by the number of columns in this case is more than two
hundred times the maximum value that can be represented by an unsigned int.
The use of @code{numel} should therefore be avoided useless it is known
it won't overflow.

Extreme care must be take with the elem method and the "()" operator,
which perform basically the same function. The reason is that if a
sparse object is non-const, then @sc{Octave} will assume that a
request for a zero element in a sparse matrix is in fact a request 
to create this element so it can be filled. Therefore a piece of
code like

@example
  SparseMatrix sm;
  @dots{}
  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      std::cerr << " (" << i << "," << j << "): " << sm(i,j) 
                << std::endl;
@end example

is a great way of turning the sparse matrix into a dense one, and a
very slow way at that since it reallocates the sparse object at each
zero element in the matrix.

An easy way of preventing the above from hapening is to create a temporary
constant version of the sparse matrix. Note that only the container for
the sparse matrix will be copied, while the actual representation of the
data will be shared between the two versions of the sparse matrix. So this
is not a costly operation. For example, the above would become

@example
  SparseMatrix sm;
  @dots{}
  const SparseMatrix tmp (sm);
  for (int j = 0; j < nc; j++)
    for (int i = 0; i < nr; i++)
      std::cerr << " (" << i << "," << j << "): " << tmp(i,j) 
                << std::endl;
@end example

Finally, as the sparse types aren't just represented as a contiguous
block of memory, the @code{fortran_vec} method of the @code{Array<T>}
is not available. It is however replaced by three seperate methods
@code{ridx}, @code{cidx} and @code{data}, that access the raw compressed
column format that the @sc{Octave} sparse matrices are stored in.
Additionally, these methods can be used in a manner similar to @code{elem},
to allow the matrix to be accessed or filled. However, in that case it is
up to the user to repect the sparse matrix compressed column format
discussed previous.

@node OctCreation, OctUse, OctDifferences, Oct-Files
@subsection Creating Spare Matrices in Oct-Files

The user has several alternatives in how to create a sparse matrix.
They can first create the data as three vectors representing the
row and column indexes and the data, and from those create the matrix.
Or alternatively, they can create a sparse matrix with the appropriate
amount of space and then fill in the values. Both techniques have their
advantages and disadvantages.

An example of how to create a small sparse matrix with the first technique
might be seen the example

@example
  int nz = 4, nr = 3, nc = 4;
  ColumnVector ridx (nz);
  ColumnVector cidx (nz);
  ColumnVector data (nz);

  ridx(0) = 0; ridx(1) = 0; ridx(2) = 1; ridx(3) = 2;
  cidx(0) = 0; cidx(1) = 1; cidx(2) = 3; cidx(3) = 3;
  data(0) = 1; data(1) = 2; data(2) = 3; data(3) = 4;

  SparseMatrix sm (data, ridx, cidx, nr, nc);
@end example

which creates the matrix given in section @ref{Storage}. Note that 
the compressed matrix format is not used at the time of the creation
of the matrix itself, however it is used internally. 

As previously mentioned, the values of the sparse matrix are stored
in increasing column-major ordering. Although the data passed by the
user does not need to respect this requirement, the pre-sorting the
data significantly speeds up the creation of the sparse matrix.

The disadvantage of this technique of creating a sparse matrix is
that there is a brief time where two copies of the data exists. Therefore
for extremely memory constrained problems this might not be the right
technique to create the sparse matrix.

The alternative is to first create the sparse matrix with the desired
number of non-zero elements and then later fill those elements in. The
easiest way to do this is 

@example 
  int nz = 4, nr = 3, nc = 4;
  SparseMatrix sm (nr, nc, nz);
  sm(0,0) = 1; sm(0,1) = 2; sm(1,3) = 3; sm(2,3) = 4;
@end example

That creates the same matrix as previously. Again, although it is not
strictly necessary, it is significantly faster if the sparse matrix is
created in this manner that the elements are added in column-major
ordering. The reason for this is that if the elements are inserted
at the end of the current list of known elements then no element
in the matrix needs to be moved to allow the new element to be
inserted. Only the column indexes need to be updated.

There are a few further points to note about this technique of creating
a sparse matrix. Firstly, it is not illegal to create a sparse matrix 
with fewer elements than are actually inserted in the matrix. Therefore

@example 
  int nz = 4, nr = 3, nc = 4;
  SparseMatrix sm (nr, nc, 0);
  sm(0,0) = 1; sm(0,1) = 2; sm(1,3) = 3; sm(2,3) = 4;
@end example

is perfectly legal. However it is a very bad idea. The reason is that 
as each new element is added to the sparse matrix the space allocated
to it is increased by reallocating the memory. This is an expensive
operation, that will significantly slow this means of creating a sparse
matrix. Furthermore, it is not illegal to create a sparse matrix with 
too much storage, so having @var{nz} above equaling 6 is also legal.
The disadvantage is that the matrix occupies more memory than strictly
needed.

It is not always easy to known the number of non-zero elements prior
to filling a matrix. For this reason the additional storage for the
sparse matrix can be removed after its creation with the
@dfn{maybe_compress} function. Furthermore, the maybe_compress can
deallocate the unused storage, but it can equally remove zero elements
from the matrix.  The removal of zero elements from the matrix is
controlled by setting the argument of the @dfn{maybe_compress} function
to be 'true'. However, the cost of removing the zeros is high because it
implies resorting the elements. Therefore, if possible it is better
is the user doesn't add the zeros in the first place. An example of
the use of @dfn{maybe_compress} is

@example
  int nz = 6, nr = 3, nc = 4;
  SparseMatrix sm1 (nr, nc, nz);
  sm1(0,0) = 1; sm1(0,1) = 2; sm1(1,3) = 3; sm1(2,3) = 4;
  sm1.maybe_compress ();  // No zero elements were added

  SparseMatrix sm2 (nr, nc, nz);
  sm2(0,0) = 1; sm2(0,1) = 2; sm(0,2) = 0; sm(1,2) = 0; 
  sm1(1,3) = 3; sm1(2,3) = 4;
  sm2.maybe_compress (true);  // Zero elements were added
@end example

The use of the @dfn{maybe_compress} function should be avoided if
possible, as it will slow the creation of the matrices.

A third means of creating a sparse matrix is to work directly with
the data in compressed row format. An example of this technique might
be

@c Note the @verbatim environment is a relatively new addition to texinfo.
@c Therefore use the @example environment and replace @, with @@, 
@c { with @{, etc

@example
  octave_value arg;
  
  @dots{}

  int nz = 6, nr = 3, nc = 4;   // Assume we know the max no nz 
  SparseMatrix sm (nr, nc, nz);
  Matrix m = arg.matrix_value ();

  int ii = 0;
  sm.cidx (0) = 0;
  for (int j = 1; j < nc; j++)
    @{
      for (int i = 0; i < nr; i++)
        @{
          double tmp = foo (m(i,j));
          if (tmp != 0.)
            @{
              sm.data(ii) = tmp;
              sm.ridx(ii) = i;
              ii++;
            @}
        @}
      sm.cidx(j+1) = ii;
   @}
  sm.maybe_mutate ();  // If don't know a-priori the final no of nz.
@end example

which is probably the most efficient means of creating the sparse matrix.

Finally, it might sometimes arise that the amount of storage initially
created is insufficient to completely store the sparse matrix. Therefore,
the method @code{change_capacity} exists to reallocate the sparse memory.
The above example would then be modified as 

@example
  octave_value arg;
  
  @dots{}

  int nz = 6, nr = 3, nc = 4;   // Assume we know the max no nz 
  SparseMatrix sm (nr, nc, nz);
  Matrix m = arg.matrix_value ();

  int ii = 0;
  sm.cidx (0) = 0;
  for (int j = 1; j < nc; j++)
    @{
      for (int i = 0; i < nr; i++)
        @{
          double tmp = foo (m(i,j));
          if (tmp != 0.)
            @{
              if (ii == nz)
                @{
                  nz += 2;   // Add 2 more elements
                  sm.change_capacity (nz);
                @}
              sm.data(ii) = tmp;
              sm.ridx(ii) = i;
              ii++;
            @}
        @}
      sm.cidx(j+1) = ii;
   @}
  sm.maybe_mutate ();  // If don't know a-priori the final no of nz.
@end example

Note that both increasing and decreasing the number of non-zero elements in
a sparse matrix is expensive, as it involves memory reallocation. Also as
parts of the matrix, though not its entirety, exist as the old and new copy
at the same time, additional memory is needed. Therefore if possible this
should be avoided.

@node OctUse, , OctCreation, Oct-Files
@subsection Using Sparse Matrices in Oct-Files

Most of the same operators and functions on sparse matrices that are
available from the @sc{Octave} are equally available with oct-files.
The basic means of extracting a sparse matrix from an @code{octave_value}
and returning them as an @code{octave_value}, can be seen in the
following example

@example
   octave_value_list retval;

   SparseMatrix sm = args(0).sparse_matrix_value ();
   SparseComplexMatrix scm = args(1).sparse_complex_matrix_value ();
   SparseBoolMatrix sbm = args(2).sparse_bool_matrix_value ();

   @dots{}

   retval(2) = sbm;
   retval(1) = scm;
   retval(0) = sm;
@end example

The conversion to an octave-value is handled by the sparse
@code{octave_value} constructors, and so no special care is needed.

@node License, Function Reference, Oct-Files, Sparse Matrices
@section Licensing of Third Party Software

There are several third party software packages used by the @sc{Octave}
sparse matrix.

@table @asis
@item COLAMD
is used for the @dfn{colamd} and @dfn{symamd} functions.

@table @asis
@item Authors
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.

@item License
Copyright @copyright{} 1998-2003 by the University of Florida.
All Rights Reserved.

THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.

Permission is hereby granted to use, copy, modify, and/or distribute
this program, provided that the Copyright, this License, and the
Availability of the original version is retained on all copies and made
accessible to the end-user of any code or package that includes COLAMD
or any modified version of COLAMD. 

@item Availability
@url{http://www.cise.ufl.edu/research/sparse/colamd/}
@end table

@item UMFPACK
is used in various places with the sparse types, including the
LU decomposition and solvers.

@table @asis
@item License
UMFPACK Version 4.4 (Jan. 28, 2005), Copyright @copyright{} 2005 by
Timothy A. Davis.  All Rights Reserved.

Your use or distribution of UMFPACK or any modified version of
UMFPACK implies that you agree to this License.

THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.

Permission is hereby granted to use or copy this program, provided
that the Copyright, this License, and the Availability of the original
version is retained on all copies.  User documentation of any code that
uses UMFPACK or any modified version of UMFPACK code must cite the
Copyright, this License, the Availability note, and 'Used by permission.'
Permission to modify the code and to distribute modified code is granted,
provided the Copyright, this License, and the Availability note are
retained, and a notice that the code was modified is included.  This
software was developed with support from the National Science Foundation,
and is provided to you free of charge.

@item Availability
@url{http://www.cise.ufl.edu/research/sparse/umfpack/}
@end table
@end table

@node Function Reference, , License, Sparse Matrices
@section Function Reference

@iftex
@subsection Functions by Category
@subsubsection Generate sparse matrix
@table @asis
@item spdiags
A generalization of the function `spdiag'.
@item speye
Returns a sparse identity matrix.
@item sprand
Generate a random sparse matrix.
@item sprandn
Generate a random sparse matrix.
@item sprandsym
@emph{Not implemented}
@end table
@subsubsection Sparse matrix conversion
@table @asis
@item full
returns a full storage matrix from a sparse one See also: sparse
@item sparse
SPARSE: create a sparse matrix
@item spconvert
This function converts for a simple sparse matrix format easily produced by other programs into Octave's internal sparse format.
@item spfind
SPFIND: a sparse version of the find operator 1.
@end table
@subsubsection Manipulate sparse matrices
@table @asis
@item issparse
Return 1 if the value of the expression EXPR is a sparse matrix.
@item nnz
returns number of non zero elements in SM See also: sparse
@item nonzeros
Returns a vector of the non-zero values of the sparse matrix S
@item nzmax
Returns the amount of storage allocated to the sparse matrix SM.
@item spalloc
Returns an empty sparse matrix of size R-by-C.
@item spfun
Compute `f(X)' for the non-zero values of X This results in a sparse matrix with the same structure as X.
@item spones
Replace the non-zero entries of X with ones.
@item spy
Plot the sparsity pattern of the sparse matrix X
@end table
@subsubsection Graph Theory
@table @asis
@item etree
Returns the elimination tree for the matrix S.
@item etreeplot
@emph{Not implemented}
@item gplot
@emph{Not implemented}
@item treelayout
@emph{Not implemented}
@item treeplot
@emph{Not implemented}
@end table
@subsubsection Sparse matrix reordering
@table @asis
@item colamd
Column approximate minimum degree permutation.
@item colperm
Returns the column permutations such that the columns of `S (:, P)' are ordered in terms of increase number of non-zero elements.
@item dmperm
Perfrom a Deulmage-Mendelsohn permutation on the sparse matrix S.
@item symamd
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.
@item symrcm
Returns the Reverse Cuthill McKee reordering of the sparse matrix S.
@end table
@subsubsection Linear algebra
@table @asis
@item cholinc
@emph{Not implemented}
@item condest
@emph{Not implemented}
@item eigs
@emph{Not implemented}
@item normest
@emph{Not implemented}
@item spdet
Compute the determinant of sparse matrix A using UMFPACK.
@item spinv
Compute the inverse of the square matrix A.
@item splu
Compute the LU decomposition of the sparse matrix A, using subroutines from UMFPACK.
@item sprank
@emph{Not implemented}
@item svds
@emph{Not implemented}
@end table
@subsubsection Iterative techniques
@table @asis
@item bicg
@emph{Not implemented}
@item bicgstab
@emph{Not implemented}
@item cgs
@emph{Not implemented}
@item gmres
@emph{Not implemented}
@item luinc
Produce the incomplete LU factorization of the sparse matrix A.
@item lsqr
@emph{Not implemented}
@item minres
@emph{Not implemented}
@item pcg
@emph{Not implemented}
@item pcr
@emph{Not implemented}
@item qmr
@emph{Not implemented}
@item symmlq
@emph{Not implemented}
@end table
@subsubsection Miscellaneous
@table @asis
@item spaugment
@emph{Not implemented}
@item spparms
Sets or displays the parameters used by the sparse solvers and factorization functions.
@item symbfact
Performs a symbolic factorization analysis on the sparse matrix S.
@item spstats
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
@item spprod
Product of elements along dimension DIM.
@item spcumprod
Cumulative product of elements along dimension DIM.
@item spcumsum
Cumulative sum of elements along dimension DIM.
@item spsum
Sum of elements along dimension DIM.
@item spsumsq
Sum of squares of elements along dimension DIM.
@item spmin
For a vector argument, return the minimum value.
@item spmax
For a vector argument, return the maximum value.
@item spatan2
Compute atan (Y / X) for corresponding sparse matrix elements of Y and X.
@item spdiag
Return a diagonal matrix with the sparse vector V on diagonal K.
@item spreshape
Return a sparse matrix with M rows and N columns whose elements are taken from the sparse matrix A.
@end table

@subsection Functions Alphabetically
@end iftex

@menu
* colamd::	Column approximate minimum degree permutation.
* colperm::	Returns the column permutations such that the columns of `S
		(:, P)' are ordered in terms of increase number of non-zero
		elements.
* dmperm::	Perfrom a Deulmage-Mendelsohn permutation on the sparse
		matrix S.
* etree::	Returns the elimination tree for the matrix S.
* full::	returns a full storage matrix from a sparse one See also:
		sparse
* issparse::	Return 1 if the value of the expression EXPR is a sparse
		matrix.
* luinc::	Produce the incomplete LU factorization of the sparse 
		A.
* nnz:: 	returns number of non zero elements in SM See also: sparse
* nonzeros::	Returns a vector of the non-zero values of the sparse
		matrix S
* nzmax::	Returns the amount of storage allocated to the sparse
		matrix SM.
* spalloc::	Returns an empty sparse matrix of size R-by-C.
* sparse::	SPARSE: create a sparse matrix
* spatan2::	Compute atan (Y / X) for corresponding sparse matrix
		elements of Y and X.
* spconvert::	This function converts for a simple sparse matrix format
		easily produced by other programs into Octave's internal
		sparse format.
* spcumprod::	Cumulative product of elements along dimension DIM.
* spcumsum::	Cumulative sum of elements along dimension DIM.
* spdet::	Compute the determinant of sparse matrix A using UMFPACK.
* spdiag::	Return a diagonal matrix with the sparse vector V on
		diagonal K.
* spdiags::	A generalization of the function `spdiag'.
* speye::	Returns a sparse identity matrix.
* spfind::	SPFIND: a sparse version of the find operator 1.
* spfun::	Compute `f(X)' for the non-zero values of X This results in
		a sparse matrix with the same structure as X.
* spinv::	Compute the inverse of the square matrix A.
* splu::	Compute the LU decomposition of the sparse matrix A, using
		subroutines from UMFPACK.
* spmax::	For a vector argument, return the maximum value.
* spmin::	For a vector argument, return the minimum value.
* spones::	Replace the non-zero entries of X with ones.
* spparms::	Sets or displays the parameters used by the sparse solvers
		and factorization functions.
* spprod::	Product of elements along dimension DIM.
* sprand::	Generate a random sparse matrix.
* sprandn::	Generate a random sparse matrix.
* spreshape::	Return a sparse matrix with M rows and N columns whose
		elements are taken from the sparse matrix A.
* spstats::	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
* spsum::	Sum of elements along dimension DIM.
* spsumsq::	Sum of squares of elements along dimension DIM.
* spy:: 	Plot the sparsity pattern of the sparse matrix X
* symamd::	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.
* symbfact::	Performs a symbolic factorization analysis on the sparse
		matrix S.
* symrcm::	Returns the Reverse Cuthill McKee reordering of the sparse
		matrix S.
@end menu

@node colamd, colperm, , Function Reference
@subsubsection colamd

@DOCSTRING(colamd)

@node colperm, dmperm, colamd, Function Reference
@subsubsection colperm

@DOCSTRING(colperm)

@node dmperm, etree, colperm, Function Reference
@subsubsection dmperm

@DOCSTRING(dmperm)

@node etree, full, dmperm, Function Reference
@subsubsection etree

@DOCSTRING(etree)

@node full, issparse, etree, Function Reference
@subsubsection full

@DOCSTRING(full)

@node issparse, luinc, full, Function Reference
@subsubsection issparse

@DOCSTRING(issparse)

@node luinc, nnz, issparse, Function Reference
@subsubsection luinc

@DOCSTRING(luinc)

@node nnz, nonzeros, luinc, Function Reference
@subsubsection nnz

@DOCSTRING(nnz)

@node nonzeros, nzmax, nnz, Function Reference
@subsubsection nonzeros

@DOCSTRING(nonzeros)

@node nzmax, spalloc, nonzeros, Function Reference
@subsubsection nzmax

@DOCSTRING(nzmax)

@node spalloc, sparse, nzmax, Function Reference
@subsubsection spalloc

@DOCSTRING(spalloc)

@node sparse, spatan2, spalloc, Function Reference
@subsubsection sparse

@DOCSTRING(sparse)

@node spatan2, spconvert, sparse, Function Reference
@subsubsection spatan2

@DOCSTRING(spatan2)

@node spconvert, spcumprod, spatan2, Function Reference
@subsubsection spconvert

@DOCSTRING(spconvert)

@node spcumprod, spcumsum, spconvert, Function Reference
@subsubsection spcumprod

@DOCSTRING(spcumprod)

@node spcumsum, spdet, spcumprod, Function Reference
@subsubsection spcumsum

@DOCSTRING(spcumsum)

@node spdet, spdiag, spcumsum, Function Reference
@subsubsection spdet

@DOCSTRING(spdet)

@node spdiag, spdiags, spdet, Function Reference
@subsubsection spdiag

@DOCSTRING(spdiag)

@node spdiags, speye, spdiag, Function Reference
@subsubsection spdiags

@DOCSTRING(spdiags)

@node speye, spfind, spdiags, Function Reference
@subsubsection speye

@DOCSTRING(speye)

@node spfind, spfun, speye, Function Reference
@subsubsection spfind

@DOCSTRING(spfind)

@node spfun, spinv, spfind, Function Reference
@subsubsection spfun

@DOCSTRING(spfun)

@node spinv, splu, spfun, Function Reference
@subsubsection spinv

@DOCSTRING(spinv)

@node splu, spmax, spinv, Function Reference
@subsubsection splu

@DOCSTRING(splu)

@node spmax, spmin, splu, Function Reference
@subsubsection spmax

@DOCSTRING(spmax)

@node spmin, spones, spmax, Function Reference
@subsubsection spmin

@DOCSTRING(spmin)

@node spones, spparms, spmin, Function Reference
@subsubsection spones

@DOCSTRING(spones)

@node spparms, spprod, spones, Function Reference
@subsubsection spparms

@DOCSTRING(spparms)

@node spprod, sprand, spparms, Function Reference
@subsubsection spprod

@DOCSTRING(spprod)

@node sprand, sprandn, spprod, Function Reference
@subsubsection sprand

@DOCSTRING(sprand)

@node sprandn, spreshape, sprand, Function Reference
@subsubsection sprandn

@DOCSTRING(sprandn)

@node spreshape, spstats, sprandn, Function Reference
@subsubsection spreshape

@DOCSTRING(spreshape)

@node spstats, spsum, spreshape, Function Reference
@subsubsection spstats

@DOCSTRING(spstats)

@node spsum, spsumsq, spstats, Function Reference
@subsubsection spsum

@DOCSTRING(spsum)

@node spsumsq, spy, spsum, Function Reference
@subsubsection spsumsq

@DOCSTRING(spsumsq)

@node spy, symamd, spsumsq, Function Reference
@subsubsection spy

@DOCSTRING(spy)

@node symamd, symbfact, spy, Function Reference
@subsubsection symamd

@DOCSTRING(symamd)

@node symbfact, symrcm, symamd, Function Reference
@subsubsection symbfact

@DOCSTRING(symbfact)

@node symrcm, , symbfact, Function Reference
@subsubsection symrcm

@DOCSTRING(symrcm)

@c Local Variables: ***
@c Mode: texinfo ***
@c End: ***