Mercurial > octave-antonio
diff doc/interpreter/sparse.txi @ 5164:57077d0ddc8e
[project @ 2005-02-25 19:55:24 by jwe]
author | jwe |
---|---|
date | Fri, 25 Feb 2005 19:55:28 +0000 |
parents | |
children | 5bdc3f24cd5f |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/doc/interpreter/sparse.txi Fri Feb 25 19:55:28 2005 +0000 @@ -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: ***