--- a/Makeconf.in
+++ b/Makeconf.in
@@ -101,13 +101,21 @@
   TMP_IF_1 = -I. -I$(srcdir)
 endif
 ifeq ($(TOPDIR),$(top_srcdir))
-  TMP_IF_2 = -I$(TOPDIR) -I$(TOPDIR)/liboctave -I$(TOPDIR)/src \
-	-I$(TOPDIR)/libcruft/misc
+  TMP_IF_2 = \
+    -I$(TOPDIR) \
+    -I$(TOPDIR)/liboctave \
+    -I$(TOPDIR)/src \
+    -I$(TOPDIR)/libcruft/misc
 else
-  TMP_IF_2 = -I$(TOPDIR) -I$(TOPDIR)/liboctave -I$(TOPDIR)/src \
-	-I$(TOPDIR)/libcruft/misc -I$(top_srcdir) \
-	-I$(top_srcdir)/liboctave -I$(top_srcdir)/src \
-	-I$(top_srcdir)/libcruft/misc
+  TMP_IF_2 = \
+    -I$(TOPDIR) \
+    -I$(TOPDIR)/liboctave \
+    -I$(TOPDIR)/src \
+    -I$(TOPDIR)/libcruft/misc \
+    -I$(top_srcdir) \
+    -I$(top_srcdir)/liboctave \
+    -I$(top_srcdir)/src \
+    -I$(top_srcdir)/libcruft/misc
 endif
 INCFLAGS = $(TMP_IF_1) $(TMP_IF_2) $(DLFCN_INCFLAGS) $(GLOB_INCFLAGS)
 

--- a/PROJECTS
+++ b/PROJECTS
@@ -35,8 +35,6 @@
 
   * Free NLP solver.
 
-  * Support for sparse matrices.
-
   * Fix CollocWt to handle Laguerre polynomials.  Make it easy to
     extend it to other polynomial types.
 
@@ -105,6 +103,104 @@
     non-square systems compatible with Matlab.  Currently, they return
     the minimum norm solution from DGELSS, which behaves differently.
 
+---------------
+Sparse Matrices:
+---------------
+
+  * Sparse Cholesky factorization for Fchol function and also for the 
+    in built polymorphic solvers.
+
+  * QR factorization functions, also for use in lssolve functions. Write
+    svds function based on this. Write sprank function based on svds.
+
+  * Once dmperm is implemented, use the technique to detect permuted
+    triangular matrices. Test the permuted triangular matrix solver code
+
+  * Accelerate the copying of the data from a sparse matrix to a banded matrix
+    in the solvers, that takes a significant portion of the computation time 
+    for narrow matrices. This is not obvious, due to the treatment of zero
+    elements in the matrix. Maybe current solution is optimal.
+
+  * Perhaps split the overly long ::solve functions up, either by the type 
+    of solver, or seperate factorization functions, so that they can be 
+    reused in each of 4 different ::solve functions.
+
+  * Sparse inverse function, based on Andy's code from octave-forge.
+
+  * Implement fourth argument to the sprand and sprandn that the leading
+    brand implements.
+
+  * Mapper functions such as real, imag, abs, etc need to be treated, either
+    with a dispatch or by changing the mapper function code.
+
+  * Write the rest of the sparse docs
+
+  * The algo in TOMS 582 is perfect for symrcm function. However, this is
+    under  the ACM license and can't be used in a GPL program.
+
+    An alternative is that PETSC is GPL compatiable and has a symrcm 
+    implemented from the original SPARSPAK. Its not clear that this is
+    legal to me as I have found no clarification of the original license
+    of SPARSPAK. As PETSC has had this code for over 10 years without
+    problem, we can perhaps assume that there is no issues. Maybe need
+    to contact PETSC people or the SPARSPAK people at uni of waterloo
+    to check issues.
+
+  * Make spalloc(r,c,n) actually create an empty sparse with n non-zero 
+    elements? This allows something like
+
+    sm = spalloc (r,c,n)
+    for j=1:c
+      for i=1:r
+        tmp = foo (i,j);
+        if (tmp != 0.)
+          sm (i,j) = tmp;
+        endif
+      endfor
+    endfor
+
+    actually make sense. Otherwise the above will cause massive amounts
+    of memory reallocation.
+
+    The fact is that this doesn't make sense in any case as the assign
+    function makes another copy of the sparse matrix. So although spalloc
+    might easily be made to have the correct behaviour, the first assign
+    will cause the matrix to be resized !!! There seems to be no simple
+    way to treat this but a complete rewrite of the sparse assignment
+    functions...
+
+  * Port the sparse testing code into the DejaGNU testing code.
+
+  * Treat the dispatching of the splu, spdet, functions, etc within octave
+    itself. This either means that classes need implementing or at the
+    minimum the octave-forge dispatch function is implemented within octave
+
+  * Other missing Functions
+      - eigs
+      - dmperm      Tim Davis is apparently working on something
+      - symmmd      Superseded by symamd
+      - colmmd      Superseded by colamd
+      - sprandsym
+      - symbfact
+      - etreeplot
+      - treeplot
+      - gplot
+      - treelayout
+      - cholinc
+      - condest
+      - normest
+      - bicg
+      - bicgstab
+      - cgs
+      - gmres
+      - lsqr
+      - minres
+      - pcg         Can this be taken from o-forge? Is it compatiable
+      - pcr         Can and should this be taken from o-forge?
+      - qmr
+      - symmlq 
+      - spaugment
+
 --------
 Graphics:
 --------

--- a/doc/ChangeLog
+++ b/doc/ChangeLog
@@ -1,3 +1,14 @@
+2005-02-25  John W. Eaton  <jwe@octave.org>
+
+	Sparse merge.
+
+	2005-01-07  David Bateman  <dbateman@free.fr>
+
+	* interpreter/sparse.txi: New file.
+	* interpreter/Makefile.in: Add it to the build.
+	* interpreter/octave.texi: Include it in the contents tables and the
+	top level document.
+
 2005-02-25  John W. Eaton  <jwe@octave.org>
 
 	* liboctave/Makefile.in (all): Also depend on Octave-FAQ.pdf.

--- a/doc/interpreter/Makefile.in
+++ b/doc/interpreter/Makefile.in
@@ -25,7 +25,7 @@
 	grammar.txi image.txi install.txi intro.txi io.txi \
 	linalg.txi matrix.txi nonlin.txi numbers.txi \
 	op-idx.txi optim.txi plot.txi poly.txi preface.txi \
-	quad.txi quaternion.txi set.txi signal.txi stats.txi \
+	quad.txi quaternion.txi set.txi signal.txi sparse.txi stats.txi \
 	stmt.txi stream.txi strings.txi struct.txi system.txi \
 	tips.txi var.txi vr-idx.txi
 

--- a/doc/interpreter/octave.texi
+++ b/doc/interpreter/octave.texi
@@ -138,6 +138,7 @@
 * Arithmetic::                  
 * Linear Algebra::              
 * Nonlinear Equations::         
+* Sparse Matrices::
 * Quadrature::                  
 * Differential Equations::      
 * Optimization::                
@@ -378,6 +379,16 @@
 * Matrix Factorizations::       
 * Functions of a Matrix::       
 
+Sparse Matrices
+
+* Basics::
+* Graph Theory::
+* Sparse Linear Algebra::
+* Iterative Techniques::
+* Oct-Files::
+* License::
+* Function Reference::
+
 Quadrature
 
 * Functions of One Variable::   
@@ -507,6 +518,7 @@
 @include arith.texi
 @include linalg.texi
 @include nonlin.texi
+@include sparse.texi
 @include quad.texi
 @include diffeq.texi
 @include optim.texi

new file mode 100644
--- /dev/null
+++ b/doc/interpreter/sparse.txi
@@ -0,0 +1,1308 @@
+@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.3 (Jan. 16, 2004), Copyright @copyright{} 2004 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 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.
+* 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, nnz, full, Function Reference
+@subsubsection issparse
+
+@DOCSTRING(issparse)
+
+@node nnz, nonzeros, issparse, 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: ***

--- a/libcruft/ChangeLog
+++ b/libcruft/ChangeLog
@@ -1,3 +1,29 @@
+2005-02-25  John W. Eaton  <jwe@octave.org>
+
+	* blas/zher.f: New file.
+
+	Sparse merge.
+
+	2005-01-13 David Bateman  <dbateman@free.fr>
+
+	* lapack/dgttrf.f lapack/dgttrs.f lapacl/zgttrf.f lapack/zgttrs.f:
+	new files
+
+	2005-01-23  David Bateman  <dbateman@free.fr>
+
+	* lapack/dgtsv.f lapack/dpbcon.f lapack/dpbtf2.f lapack/dpbtrf.f 
+	lapack/dpbtrs.f lapack/dptsv.f lapack/dpttrf.f lapack/dpttrs.f 
+	lapack/dptts2.f lapack/zgtsv.f lapack/zpbcon.f lapack/zpbtf2.f 
+	lapack/zpbtrf.f lapack/zpbtrs.f lapack/zptsv.f lapack/zpttrf.f 
+	lapack/zpttrs.f lapck/zptts2.f: New files.
+
+	2004-12-29  John W. Eaton  <jwe@octave.org>
+
+	* blas/zbtsv.f: New file.
+	* lapack/dgbcon.f, lapack/dgbtrf.f, lapack/dgbtrs.f,
+	lapack/dlatbs.f, lapack/zgbcon.f, lapack/zgbtf2.f,
+	lapack/zgbtrf.f, lapack/zgbtrs.f, lapack/zlatbs.f: New files.
+
 2005-02-10  John W. Eaton  <jwe@octave.org>
 
 	* misc/cquit.c (octave_signal_caught): New global variable.

new file mode 100644
--- /dev/null
+++ b/libcruft/blas/zher.f
@@ -0,0 +1,212 @@
+      SUBROUTINE ZHER  ( UPLO, N, ALPHA, X, INCX, A, LDA )
+*     .. Scalar Arguments ..
+      DOUBLE PRECISION   ALPHA
+      INTEGER            INCX, LDA, N
+      CHARACTER*1        UPLO
+*     .. Array Arguments ..
+      COMPLEX*16         A( LDA, * ), X( * )
+*     ..
+*
+*  Purpose
+*  =======
+*
+*  ZHER   performs the hermitian rank 1 operation
+*
+*     A := alpha*x*conjg( x' ) + A,
+*
+*  where alpha is a real scalar, x is an n element vector and A is an
+*  n by n hermitian matrix.
+*
+*  Parameters
+*  ==========
+*
+*  UPLO   - CHARACTER*1.
+*           On entry, UPLO specifies whether the upper or lower
+*           triangular part of the array A is to be referenced as
+*           follows:
+*
+*              UPLO = 'U' or 'u'   Only the upper triangular part of A
+*                                  is to be referenced.
+*
+*              UPLO = 'L' or 'l'   Only the lower triangular part of A
+*                                  is to be referenced.
+*
+*           Unchanged on exit.
+*
+*  N      - INTEGER.
+*           On entry, N specifies the order of the matrix A.
+*           N must be at least zero.
+*           Unchanged on exit.
+*
+*  ALPHA  - DOUBLE PRECISION.
+*           On entry, ALPHA specifies the scalar alpha.
+*           Unchanged on exit.
+*
+*  X      - COMPLEX*16       array of dimension at least
+*           ( 1 + ( n - 1 )*abs( INCX ) ).
+*           Before entry, the incremented array X must contain the n
+*           element vector x.
+*           Unchanged on exit.
+*
+*  INCX   - INTEGER.
+*           On entry, INCX specifies the increment for the elements of
+*           X. INCX must not be zero.
+*           Unchanged on exit.
+*
+*  A      - COMPLEX*16       array of DIMENSION ( LDA, n ).
+*           Before entry with  UPLO = 'U' or 'u', the leading n by n
+*           upper triangular part of the array A must contain the upper
+*           triangular part of the hermitian matrix and the strictly
+*           lower triangular part of A is not referenced. On exit, the
+*           upper triangular part of the array A is overwritten by the
+*           upper triangular part of the updated matrix.
+*           Before entry with UPLO = 'L' or 'l', the leading n by n
+*           lower triangular part of the array A must contain the lower
+*           triangular part of the hermitian matrix and the strictly
+*           upper triangular part of A is not referenced. On exit, the
+*           lower triangular part of the array A is overwritten by the
+*           lower triangular part of the updated matrix.
+*           Note that the imaginary parts of the diagonal elements need
+*           not be set, they are assumed to be zero, and on exit they
+*           are set to zero.
+*
+*  LDA    - INTEGER.
+*           On entry, LDA specifies the first dimension of A as declared
+*           in the calling (sub) program. LDA must be at least
+*           max( 1, n ).
+*           Unchanged on exit.
+*
+*
+*  Level 2 Blas routine.
+*
+*  -- Written on 22-October-1986.
+*     Jack Dongarra, Argonne National Lab.
+*     Jeremy Du Croz, Nag Central Office.
+*     Sven Hammarling, Nag Central Office.
+*     Richard Hanson, Sandia National Labs.
+*
+*
+*     .. Parameters ..
+      COMPLEX*16         ZERO
+      PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
+*     .. Local Scalars ..
+      COMPLEX*16         TEMP
+      INTEGER            I, INFO, IX, J, JX, KX
+*     .. External Functions ..
+      LOGICAL            LSAME
+      EXTERNAL           LSAME
+*     .. External Subroutines ..
+      EXTERNAL           XERBLA
+*     .. Intrinsic Functions ..
+      INTRINSIC          DCONJG, MAX, DBLE
+*     ..
+*     .. Executable Statements ..
+*
+*     Test the input parameters.
+*
+      INFO = 0
+      IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
+     $         .NOT.LSAME( UPLO, 'L' )      )THEN
+         INFO = 1
+      ELSE IF( N.LT.0 )THEN
+         INFO = 2
+      ELSE IF( INCX.EQ.0 )THEN
+         INFO = 5
+      ELSE IF( LDA.LT.MAX( 1, N ) )THEN
+         INFO = 7
+      END IF
+      IF( INFO.NE.0 )THEN
+         CALL XERBLA( 'ZHER  ', INFO )
+         RETURN
+      END IF
+*
+*     Quick return if possible.
+*
+      IF( ( N.EQ.0 ).OR.( ALPHA.EQ.DBLE( ZERO ) ) )
+     $   RETURN
+*
+*     Set the start point in X if the increment is not unity.
+*
+      IF( INCX.LE.0 )THEN
+         KX = 1 - ( N - 1 )*INCX
+      ELSE IF( INCX.NE.1 )THEN
+         KX = 1
+      END IF
+*
+*     Start the operations. In this version the elements of A are
+*     accessed sequentially with one pass through the triangular part
+*     of A.
+*
+      IF( LSAME( UPLO, 'U' ) )THEN
+*
+*        Form  A  when A is stored in upper triangle.
+*
+         IF( INCX.EQ.1 )THEN
+            DO 20, J = 1, N
+               IF( X( J ).NE.ZERO )THEN
+                  TEMP = ALPHA*DCONJG( X( J ) )
+                  DO 10, I = 1, J - 1
+                     A( I, J ) = A( I, J ) + X( I )*TEMP
+   10             CONTINUE
+                  A( J, J ) = DBLE( A( J, J ) ) + DBLE( X( J )*TEMP )
+               ELSE
+                  A( J, J ) = DBLE( A( J, J ) )
+               END IF
+   20       CONTINUE
+         ELSE
+            JX = KX
+            DO 40, J = 1, N
+               IF( X( JX ).NE.ZERO )THEN
+                  TEMP = ALPHA*DCONJG( X( JX ) )
+                  IX   = KX
+                  DO 30, I = 1, J - 1
+                     A( I, J ) = A( I, J ) + X( IX )*TEMP
+                     IX        = IX        + INCX
+   30             CONTINUE
+                  A( J, J ) = DBLE( A( J, J ) ) + DBLE( X( JX )*TEMP )
+               ELSE
+                  A( J, J ) = DBLE( A( J, J ) )
+               END IF
+               JX = JX + INCX
+   40       CONTINUE
+         END IF
+      ELSE
+*
+*        Form  A  when A is stored in lower triangle.
+*
+         IF( INCX.EQ.1 )THEN
+            DO 60, J = 1, N
+               IF( X( J ).NE.ZERO )THEN
+                  TEMP      = ALPHA*DCONJG( X( J ) )
+                  A( J, J ) = DBLE( A( J, J ) ) + DBLE( TEMP*X( J ) )
+                  DO 50, I = J + 1, N
+                     A( I, J ) = A( I, J ) + X( I )*TEMP
+   50             CONTINUE
+               ELSE
+                  A( J, J ) = DBLE( A( J, J ) )
+               END IF
+   60       CONTINUE
+         ELSE
+            JX = KX
+            DO 80, J = 1, N
+               IF( X( JX ).NE.ZERO )THEN
+                  TEMP      = ALPHA*DCONJG( X( JX ) )
+                  A( J, J ) = DBLE( A( J, J ) ) + DBLE( TEMP*X( JX ) )
+                  IX        = JX
+                  DO 70, I = J + 1, N
+                     IX        = IX        + INCX
+                     A( I, J ) = A( I, J ) + X( IX )*TEMP
+   70             CONTINUE
+               ELSE
+                  A( J, J ) = DBLE( A( J, J ) )
+               END IF
+               JX = JX + INCX
+   80       CONTINUE
+         END IF
+      END IF
+*
+      RETURN
+*
+*     End of ZHER  .
+*
+      END