Mercurial > hg > octave-nkf
annotate doc/interpreter/sparse.txi @ 9088:77e71f3da3d6
Fix documentation image printing under new development code
Printed images are now clipped to image size rather than to papersize of 8.5x11
Images are scaled to 4 inches and centered in display
pdf printing is sub-optimal as texi2pdf picks the png rendering (bitmap) rather
than the pdf rendering (vector).
| author | Rik <rdrider0-list@yahoo.com> |
|---|---|
| date | Sat, 04 Apr 2009 14:28:22 -0700 |
| parents | be150a172010 |
| children | 923c7cb7f13f |
| rev | line source |
|---|---|
| 8920 | 1 @c Copyright (C) 2004, 2005, 2006, 2007, 2008, 2009 David Bateman |
| 7018 | 2 @c |
| 3 @c This file is part of Octave. | |
| 4 @c | |
| 5 @c Octave is free software; you can redistribute it and/or modify it | |
| 6 @c under the terms of the GNU General Public License as published by the | |
| 7 @c Free Software Foundation; either version 3 of the License, or (at | |
| 8 @c your option) any later version. | |
| 9 @c | |
| 10 @c Octave is distributed in the hope that it will be useful, but WITHOUT | |
| 11 @c ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
| 12 @c FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
| 13 @c for more details. | |
| 14 @c | |
| 15 @c You should have received a copy of the GNU General Public License | |
| 16 @c along with Octave; see the file COPYING. If not, see | |
| 17 @c <http://www.gnu.org/licenses/>. | |
| 5164 | 18 |
| 5648 | 19 @ifhtml |
| 20 @set htmltex | |
| 21 @end ifhtml | |
| 22 @iftex | |
| 23 @set htmltex | |
| 24 @end iftex | |
| 25 | |
| 5164 | 26 @node Sparse Matrices |
| 27 @chapter Sparse Matrices | |
| 28 | |
| 29 @menu | |
|
8817
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
30 * Basics:: Creation and Manipulation of Sparse Matrices |
|
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
31 * Sparse Linear Algebra:: Linear Algebra on Sparse Matrices |
|
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
32 * Iterative Techniques:: Iterative Techniques |
|
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
33 * Real Life Example:: Using Sparse Matrices |
| 5164 | 34 @end menu |
| 35 | |
|
8817
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
36 @node Basics |
| 5164 | 37 @section The Creation and Manipulation of Sparse Matrices |
| 38 | |
| 39 The size of mathematical problems that can be treated at any particular | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
40 time is generally limited by the available computing resources. Both, |
| 5164 | 41 the speed of the computer and its available memory place limitation on |
| 42 the problem size. | |
| 43 | |
| 5506 | 44 There are many classes of mathematical problems which give rise to |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
45 matrices, where a large number of the elements are zero. In this case |
| 5164 | 46 it makes sense to have a special matrix type to handle this class of |
| 47 problems where only the non-zero elements of the matrix are | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
48 stored. Not only does this reduce the amount of memory to store the |
| 5164 | 49 matrix, but it also means that operations on this type of matrix can |
| 50 take advantage of the a-priori knowledge of the positions of the | |
| 51 non-zero elements to accelerate their calculations. | |
| 52 | |
| 53 A matrix type that stores only the non-zero elements is generally called | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
54 sparse. It is the purpose of this document to discuss the basics of the |
| 5164 | 55 storage and creation of sparse matrices and the fundamental operations |
| 56 on them. | |
| 57 | |
| 58 @menu | |
|
8817
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
59 * Storage of Sparse Matrices:: |
|
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
60 * Creating Sparse Matrices:: |
|
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
61 * Information:: |
|
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
62 * Operators and Functions:: |
| 5164 | 63 @end menu |
| 64 | |
|
8817
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
65 @node Storage of Sparse Matrices |
| 5164 | 66 @subsection Storage of Sparse Matrices |
| 67 | |
| 68 It is not strictly speaking necessary for the user to understand how | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
69 sparse matrices are stored. However, such an understanding will help |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
70 to get an understanding of the size of sparse matrices. Understanding |
| 5164 | 71 the storage technique is also necessary for those users wishing to |
| 72 create their own oct-files. | |
| 73 | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
74 There are many different means of storing sparse matrix data. What all |
| 5648 | 75 of the methods have in common is that they attempt to reduce the complexity |
| 5164 | 76 and storage given a-priori knowledge of the particular class of problems |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
77 that will be solved. A good summary of the available techniques for storing |
| 5164 | 78 sparse matrix is given by Saad @footnote{Youcef Saad "SPARSKIT: A basic toolkit |
| 79 for sparse matrix computation", 1994, | |
| 6620 | 80 @url{http://www-users.cs.umn.edu/~saad/software/SPARSKIT/paper.ps}}. |
| 5164 | 81 With full matrices, knowledge of the point of an element of the matrix |
| 82 within the matrix is implied by its position in the computers memory. | |
| 83 However, this is not the case for sparse matrices, and so the positions | |
| 84 of the non-zero elements of the matrix must equally be stored. | |
| 85 | |
| 86 An obvious way to do this is by storing the elements of the matrix as | |
| 5506 | 87 triplets, with two elements being their position in the array |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
88 (rows and column) and the third being the data itself. This is conceptually |
| 5164 | 89 easy to grasp, but requires more storage than is strictly needed. |
| 90 | |
| 5648 | 91 The storage technique used within Octave is the compressed column |
| 5164 | 92 format. In this format the position of each element in a row and the |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
93 data are stored as previously. However, if we assume that all elements |
| 5164 | 94 in the same column are stored adjacent in the computers memory, then |
| 95 we only need to store information on the number of non-zero elements | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
96 in each column, rather than their positions. Thus assuming that the |
| 5164 | 97 matrix has more non-zero elements than there are columns in the |
| 98 matrix, we win in terms of the amount of memory used. | |
| 99 | |
| 100 In fact, the column index contains one more element than the number of | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
101 columns, with the first element always being zero. The advantage of |
| 7001 | 102 this is a simplification in the code, in that there is no special case |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
103 for the first or last columns. A short example, demonstrating this in |
| 5164 | 104 C is. |
| 105 | |
| 106 @example | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
107 @group |
| 5164 | 108 for (j = 0; j < nc; j++) |
| 109 for (i = cidx (j); i < cidx(j+1); i++) | |
| 5648 | 110 printf ("non-zero element (%i,%i) is %d\n", |
| 111 ridx(i), j, data(i)); | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
112 @end group |
| 5164 | 113 @end example |
| 114 | |
| 115 A clear understanding might be had by considering an example of how the | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
116 above applies to an example matrix. Consider the matrix |
| 5164 | 117 |
| 118 @example | |
| 119 @group | |
| 120 1 2 0 0 | |
| 121 0 0 0 3 | |
| 122 0 0 0 4 | |
| 123 @end group | |
| 124 @end example | |
| 125 | |
| 126 The non-zero elements of this matrix are | |
| 127 | |
| 128 @example | |
| 129 @group | |
| 130 (1, 1) @result{} 1 | |
| 131 (1, 2) @result{} 2 | |
| 132 (2, 4) @result{} 3 | |
| 133 (3, 4) @result{} 4 | |
| 134 @end group | |
| 135 @end example | |
| 136 | |
| 137 This will be stored as three vectors @var{cidx}, @var{ridx} and @var{data}, | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
138 representing the column indexing, row indexing and data respectively. The |
| 5164 | 139 contents of these three vectors for the above matrix will be |
| 140 | |
| 141 @example | |
| 142 @group | |
| 5506 | 143 @var{cidx} = [0, 1, 2, 2, 4] |
| 5164 | 144 @var{ridx} = [0, 0, 1, 2] |
| 145 @var{data} = [1, 2, 3, 4] | |
| 146 @end group | |
| 147 @end example | |
| 148 | |
| 149 Note that this is the representation of these elements with the first row | |
| 5648 | 150 and column assumed to start at zero, while in Octave itself the row and |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
151 column indexing starts at one. Thus the number of elements in the |
| 5164 | 152 @var{i}-th column is given by @code{@var{cidx} (@var{i} + 1) - |
| 153 @var{cidx} (@var{i})}. | |
| 154 | |
| 5648 | 155 Although Octave uses a compressed column format, it should be noted |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
156 that compressed row formats are equally possible. However, in the |
| 5648 | 157 context of mixed operations between mixed sparse and dense matrices, |
| 158 it makes sense that the elements of the sparse matrices are in the | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
159 same order as the dense matrices. Octave stores dense matrices in |
| 5648 | 160 column major ordering, and so sparse matrices are equally stored in |
| 161 this manner. | |
| 5164 | 162 |
| 5324 | 163 A further constraint on the sparse matrix storage used by Octave is that |
| 5164 | 164 all elements in the rows are stored in increasing order of their row |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
165 index, which makes certain operations faster. However, it imposes |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
166 the need to sort the elements on the creation of sparse matrices. Having |
| 7001 | 167 disordered elements is potentially an advantage in that it makes operations |
| 5164 | 168 such as concatenating two sparse matrices together easier and faster, however |
| 169 it adds complexity and speed problems elsewhere. | |
| 170 | |
|
8817
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
171 @node Creating Sparse Matrices |
| 5164 | 172 @subsection Creating Sparse Matrices |
| 173 | |
| 174 There are several means to create sparse matrix. | |
| 175 | |
| 176 @table @asis | |
| 177 @item Returned from a function | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
178 There are many functions that directly return sparse matrices. These include |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7505
diff
changeset
|
179 @dfn{speye}, @dfn{sprand}, @dfn{diag}, etc. |
| 5164 | 180 @item Constructed from matrices or vectors |
| 181 The function @dfn{sparse} allows a sparse matrix to be constructed from | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
182 three vectors representing the row, column and data. Alternatively, the |
| 5164 | 183 function @dfn{spconvert} uses a three column matrix format to allow easy |
| 184 importation of data from elsewhere. | |
| 185 @item Created and then filled | |
| 186 The function @dfn{sparse} or @dfn{spalloc} can be used to create an empty | |
| 187 matrix that is then filled by the user | |
| 188 @item From a user binary program | |
| 189 The user can directly create the sparse matrix within an oct-file. | |
| 190 @end table | |
| 191 | |
| 192 There are several basic functions to return specific sparse | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
193 matrices. For example the sparse identity matrix, is a matrix that is |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
194 often needed. It therefore has its own function to create it as |
| 5164 | 195 @code{speye (@var{n})} or @code{speye (@var{r}, @var{c})}, which |
| 196 creates an @var{n}-by-@var{n} or @var{r}-by-@var{c} sparse identity | |
| 197 matrix. | |
| 198 | |
| 199 Another typical sparse matrix that is often needed is a random distribution | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
200 of random elements. The functions @dfn{sprand} and @dfn{sprandn} perform |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
201 this for uniform and normal random distributions of elements. They have exactly |
| 5506 | 202 the same calling convention, where @code{sprand (@var{r}, @var{c}, @var{d})}, |
| 203 creates an @var{r}-by-@var{c} sparse matrix with a density of filled | |
| 5164 | 204 elements of @var{d}. |
| 205 | |
| 7001 | 206 Other functions of interest that directly create sparse matrices, are |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7505
diff
changeset
|
207 @dfn{diag} or its generalization @dfn{spdiags}, that can take the |
| 5164 | 208 definition of the diagonals of the matrix and create the sparse matrix |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
209 that corresponds to this. For example |
| 5164 | 210 |
| 211 @example | |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7505
diff
changeset
|
212 s = diag (sparse(randn(1,n)), -1); |
| 5164 | 213 @end example |
| 214 | |
| 215 creates a sparse (@var{n}+1)-by-(@var{n}+1) sparse matrix with a single | |
| 216 diagonal defined. | |
| 217 | |
| 6620 | 218 |
| 219 @DOCSTRING(spdiags) | |
| 220 | |
| 221 @DOCSTRING(speye) | |
| 222 | |
| 223 @DOCSTRING(spfun) | |
| 224 | |
| 225 @DOCSTRING(spmax) | |
| 226 | |
| 227 @DOCSTRING(spmin) | |
| 228 | |
| 229 @DOCSTRING(spones) | |
| 230 | |
| 231 @DOCSTRING(sprand) | |
| 232 | |
| 233 @DOCSTRING(sprandn) | |
| 234 | |
| 235 @DOCSTRING(sprandsym) | |
| 236 | |
| 5164 | 237 The recommended way for the user to create a sparse matrix, is to create |
| 5648 | 238 two vectors containing the row and column index of the data and a third |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
239 vector of the same size containing the data to be stored. For example |
| 5164 | 240 |
| 241 @example | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
242 @group |
| 6421 | 243 ri = ci = d = []; |
| 244 for j = 1:c | |
| 245 ri = [ri; randperm(r)(1:n)']; | |
| 246 ci = [ci; j*ones(n,1)]; | |
| 247 d = [d; rand(n,1)]; | |
| 248 endfor | |
| 249 s = sparse (ri, ci, d, r, c); | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
250 @end group |
| 5164 | 251 @end example |
| 252 | |
| 6421 | 253 creates an @var{r}-by-@var{c} sparse matrix with a random distribution |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
254 of @var{n} (<@var{r}) elements per column. The elements of the vectors |
| 6421 | 255 do not need to be sorted in any particular order as Octave will sort |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
256 them prior to storing the data. However, pre-sorting the data will |
| 6421 | 257 make the creation of the sparse matrix faster. |
| 5164 | 258 |
| 259 The function @dfn{spconvert} takes a three or four column real matrix. | |
| 260 The first two columns represent the row and column index respectively and | |
| 261 the third and four columns, the real and imaginary parts of the sparse | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
262 matrix. The matrix can contain zero elements and the elements can be |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
263 sorted in any order. Adding zero elements is a convenient way to define |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
264 the size of the sparse matrix. For example |
| 5164 | 265 |
| 266 @example | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
267 @group |
| 5164 | 268 s = spconvert ([1 2 3 4; 1 3 4 4; 1 2 3 0]') |
| 269 @result{} Compressed Column Sparse (rows=4, cols=4, nnz=3) | |
| 270 (1 , 1) -> 1 | |
| 271 (2 , 3) -> 2 | |
| 272 (3 , 4) -> 3 | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
273 @end group |
| 5164 | 274 @end example |
| 275 | |
| 276 An example of creating and filling a matrix might be | |
| 277 | |
| 278 @example | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
279 @group |
| 5164 | 280 k = 5; |
| 281 nz = r * k; | |
| 282 s = spalloc (r, c, nz) | |
| 283 for j = 1:c | |
| 284 idx = randperm (r); | |
| 5648 | 285 s (:, j) = [zeros(r - k, 1); ... |
| 286 rand(k, 1)] (idx); | |
| 5164 | 287 endfor |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
288 @end group |
| 5164 | 289 @end example |
| 290 | |
| 5324 | 291 It should be noted, that due to the way that the Octave |
| 5164 | 292 assignment functions are written that the assignment will reallocate |
| 5506 | 293 the memory used by the sparse matrix at each iteration of the above loop. |
| 294 Therefore the @dfn{spalloc} function ignores the @var{nz} argument and | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
295 does not preassign the memory for the matrix. Therefore, it is vitally |
| 5648 | 296 important that code using to above structure should be vectorized |
| 297 as much as possible to minimize the number of assignments and reduce the | |
| 5164 | 298 number of memory allocations. |
| 299 | |
| 6620 | 300 @DOCSTRING(full) |
| 301 | |
| 302 @DOCSTRING(spalloc) | |
| 303 | |
| 304 @DOCSTRING(sparse) | |
| 305 | |
| 306 @DOCSTRING(spconvert) | |
| 307 | |
|
8106
8a42498edb30
Clarify doc for sparse function
David Bateman <dbateman@free.fr>
parents:
7681
diff
changeset
|
308 The above problem of memory reallocation can be avoided in |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
309 oct-files. However, the construction of a sparse matrix from an oct-file |
| 8828 | 310 is more complex than can be discussed here, and |
|
8106
8a42498edb30
Clarify doc for sparse function
David Bateman <dbateman@free.fr>
parents:
7681
diff
changeset
|
311 you are referred to chapter @ref{Dynamically Linked Functions}, to have |
|
8a42498edb30
Clarify doc for sparse function
David Bateman <dbateman@free.fr>
parents:
7681
diff
changeset
|
312 a full description of the techniques involved. |
| 5164 | 313 |
|
8817
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
314 @node Information |
| 5648 | 315 @subsection Finding out Information about Sparse Matrices |
| 316 | |
| 317 There are a number of functions that allow information concerning | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
318 sparse matrices to be obtained. The most basic of these is |
| 5648 | 319 @dfn{issparse} that identifies whether a particular Octave object is |
| 320 in fact a sparse matrix. | |
| 321 | |
| 322 Another very basic function is @dfn{nnz} that returns the number of | |
| 323 non-zero entries there are in a sparse matrix, while the function | |
| 324 @dfn{nzmax} returns the amount of storage allocated to the sparse | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
325 matrix. Note that Octave tends to crop unused memory at the first |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
326 opportunity for sparse objects. There are some cases of user created |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
327 sparse objects where the value returned by @dfn{nzmax} will not be |
| 5648 | 328 the same as @dfn{nnz}, but in general they will give the same |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
329 result. The function @dfn{spstats} returns some basic statistics on |
| 5648 | 330 the columns of a sparse matrix including the number of elements, the |
| 331 mean and the variance of each column. | |
| 332 | |
| 6620 | 333 @DOCSTRING(issparse) |
| 334 | |
| 335 @DOCSTRING(nnz) | |
| 336 | |
| 337 @DOCSTRING(nonzeros) | |
| 338 | |
| 339 @DOCSTRING(nzmax) | |
| 340 | |
| 341 @DOCSTRING(spstats) | |
| 342 | |
| 5648 | 343 When solving linear equations involving sparse matrices Octave |
| 344 determines the means to solve the equation based on the type of the | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
345 matrix as discussed in @ref{Sparse Linear Algebra}. Octave probes the |
| 5648 | 346 matrix type when the div (/) or ldiv (\) operator is first used with |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
347 the matrix and then caches the type. However the @dfn{matrix_type} |
| 5648 | 348 function can be used to determine the type of the sparse matrix prior |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
349 to use of the div or ldiv operators. For example |
| 5648 | 350 |
| 351 @example | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
352 @group |
| 5648 | 353 a = tril (sprandn(1024, 1024, 0.02), -1) ... |
| 354 + speye(1024); | |
| 355 matrix_type (a); | |
| 356 ans = Lower | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
357 @end group |
| 5648 | 358 @end example |
| 359 | |
| 360 show that Octave correctly determines the matrix type for lower | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
361 triangular matrices. @dfn{matrix_type} can also be used to force |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
362 the type of a matrix to be a particular type. For example |
| 5648 | 363 |
| 364 @example | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
365 @group |
| 5648 | 366 a = matrix_type (tril (sprandn (1024, ... |
| 367 1024, 0.02), -1) + speye(1024), 'Lower'); | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
368 @end group |
| 5648 | 369 @end example |
| 370 | |
| 371 This allows the cost of determining the matrix type to be | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
372 avoided. However, incorrectly defining the matrix type will result in |
| 5648 | 373 incorrect results from solutions of linear equations, and so it is |
| 374 entirely the responsibility of the user to correctly identify the | |
| 375 matrix type | |
| 376 | |
| 377 There are several graphical means of finding out information about | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
378 sparse matrices. The first is the @dfn{spy} command, which displays |
| 5648 | 379 the structure of the non-zero elements of the |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
380 matrix. @xref{fig:spmatrix}, for an example of the use of |
| 5704 | 381 @dfn{spy}. More advanced graphical information can be obtained with the |
| 5648 | 382 @dfn{treeplot}, @dfn{etreeplot} and @dfn{gplot} commands. |
| 383 | |
| 384 @float Figure,fig:spmatrix | |
|
9088
77e71f3da3d6
Fix documentation image printing under new development code
Rik <rdrider0-list@yahoo.com>
parents:
9066
diff
changeset
|
385 @center @image{spmatrix,4in} |
| 5648 | 386 @caption{Structure of simple sparse matrix.} |
| 387 @end float | |
| 388 | |
| 389 One use of sparse matrices is in graph theory, where the | |
| 7001 | 390 interconnections between nodes are represented as an adjacency |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
391 matrix. That is, if the i-th node in a graph is connected to the j-th |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
392 node. Then the ij-th node (and in the case of undirected graphs the |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
393 ji-th node) of the sparse adjacency matrix is non-zero. If each node |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
394 is then associated with a set of coordinates, then the @dfn{gplot} |
| 5648 | 395 command can be used to graphically display the interconnections |
| 396 between nodes. | |
| 397 | |
| 398 As a trivial example of the use of @dfn{gplot}, consider the example | |
| 399 | |
| 400 @example | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
401 @group |
| 5648 | 402 A = sparse([2,6,1,3,2,4,3,5,4,6,1,5], |
| 403 [1,1,2,2,3,3,4,4,5,5,6,6],1,6,6); | |
| 404 xy = [0,4,8,6,4,2;5,0,5,7,5,7]'; | |
| 405 gplot(A,xy) | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
406 @end group |
| 5648 | 407 @end example |
| 408 | |
| 409 which creates an adjacency matrix @code{A} where node 1 is connected | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
410 to nodes 2 and 6, node 2 with nodes 1 and 3, etc. The coordinates of |
| 5648 | 411 the nodes are given in the n-by-2 matrix @code{xy}. |
| 412 @ifset htmltex | |
| 413 @xref{fig:gplot}. | |
| 414 | |
| 415 @float Figure,fig:gplot | |
|
9088
77e71f3da3d6
Fix documentation image printing under new development code
Rik <rdrider0-list@yahoo.com>
parents:
9066
diff
changeset
|
416 @center @image{gplot,4in} |
| 5648 | 417 @caption{Simple use of the @dfn{gplot} command.} |
| 418 @end float | |
| 419 @end ifset | |
| 420 | |
| 421 The dependencies between the nodes of a Cholesky factorization can be | |
| 422 calculated in linear time without explicitly needing to calculate the | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
423 Cholesky factorization by the @code{etree} command. This command |
| 5648 | 424 returns the elimination tree of the matrix and can be displayed |
| 425 graphically by the command @code{treeplot(etree(A))} if @code{A} is | |
| 426 symmetric or @code{treeplot(etree(A+A'))} otherwise. | |
| 427 | |
| 6620 | 428 @DOCSTRING(spy) |
| 429 | |
| 430 @DOCSTRING(etree) | |
| 431 | |
| 432 @DOCSTRING(etreeplot) | |
| 433 | |
| 434 @DOCSTRING(gplot) | |
| 435 | |
| 436 @DOCSTRING(treeplot) | |
| 437 | |
|
8817
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
438 @DOCSTRING(treelayout) |
|
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
439 |
|
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
440 @node Operators and Functions |
| 5164 | 441 @subsection Basic Operators and Functions on Sparse Matrices |
| 442 | |
| 443 @menu | |
|
8817
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
444 * Sparse Functions:: |
|
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
445 * Return Types of Operators and Functions:: |
|
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
446 * Mathematical Considerations:: |
| 5164 | 447 @end menu |
| 448 | |
|
8817
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
449 @node Sparse Functions |
| 5648 | 450 @subsubsection Sparse Functions |
| 451 | |
| 452 An important consideration in the use of the sparse functions of | |
| 453 Octave is that many of the internal functions of Octave, such as | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
454 @dfn{diag}, cannot accept sparse matrices as an input. The sparse |
| 5648 | 455 implementation in Octave therefore uses the @dfn{dispatch} |
| 456 function to overload the normal Octave functions with equivalent | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
457 functions that work with sparse matrices. However, at any time the |
| 5648 | 458 sparse matrix specific version of the function can be used by |
| 459 explicitly calling its function name. | |
| 460 | |
| 6620 | 461 The table below lists all of the sparse functions of Octave. Note that |
| 7001 | 462 the names of the |
| 463 specific sparse forms of the functions are typically the same as | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
464 the general versions with a @dfn{sp} prefix. In the table below, and the |
| 6620 | 465 rest of this article the specific sparse versions of the functions are |
| 466 used. | |
| 467 | |
| 468 @c Table includes in comments the missing sparse functions | |
| 5648 | 469 |
| 470 @table @asis | |
| 471 @item Generate sparse matrices: | |
| 472 @dfn{spalloc}, @dfn{spdiags}, @dfn{speye}, @dfn{sprand}, | |
| 473 @dfn{sprandn}, @dfn{sprandsym} | |
| 474 | |
| 475 @item Sparse matrix conversion: | |
|
7505
f5005d9510f4
Remove dispatched sparse functions and treat in the generic versions of the functions
David Bateman <dbateman@free.fr>
parents:
7455
diff
changeset
|
476 @dfn{full}, @dfn{sparse}, @dfn{spconvert} |
| 5648 | 477 |
| 478 @item Manipulate sparse matrices | |
| 479 @dfn{issparse}, @dfn{nnz}, @dfn{nonzeros}, @dfn{nzmax}, | |
| 6620 | 480 @dfn{spfun}, @dfn{spones}, @dfn{spy} |
| 5164 | 481 |
| 5648 | 482 @item Graph Theory: |
| 483 @dfn{etree}, @dfn{etreeplot}, @dfn{gplot}, | |
| 6620 | 484 @dfn{treeplot} |
| 485 @c @dfn{treelayout} | |
| 5648 | 486 |
| 487 @item Sparse matrix reordering: | |
| 7619 | 488 @dfn{amd}, @dfn{ccolamd}, @dfn{colamd}, @dfn{colperm}, @dfn{csymamd}, |
| 6620 | 489 @dfn{dmperm}, @dfn{symamd}, @dfn{randperm}, @dfn{symrcm} |
| 5648 | 490 |
| 491 @item Linear algebra: | |
|
8417
654bcfb937bf
Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
8347
diff
changeset
|
492 @dfn{condest}, @dfn{eigs}, @dfn{matrix_type}, @dfn{normest}, @dfn{sprank}, |
|
654bcfb937bf
Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
8347
diff
changeset
|
493 @dfn{spaugment}, @dfn{svds} |
| 5648 | 494 |
| 495 @item Iterative techniques: | |
| 6620 | 496 @dfn{luinc}, @dfn{pcg}, @dfn{pcr} |
| 497 @c @dfn{bicg}, @dfn{bicgstab}, @dfn{cholinc}, @dfn{cgs}, @dfn{gmres}, | |
| 498 @c @dfn{lsqr}, @dfn{minres}, @dfn{qmr}, @dfn{symmlq} | |
| 5648 | 499 |
| 500 @item Miscellaneous: | |
|
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7505
diff
changeset
|
501 @dfn{spparms}, @dfn{symbfact}, @dfn{spstats} |
| 5648 | 502 @end table |
| 503 | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
504 In addition all of the standard Octave mapper functions (i.e., basic |
|
9088
77e71f3da3d6
Fix documentation image printing under new development code
Rik <rdrider0-list@yahoo.com>
parents:
9066
diff
changeset
|
505 math functions that take a single argument) such as @dfn{abs}, etc. |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
506 can accept sparse matrices. The reader is referred to the documentation |
| 5648 | 507 supplied with these functions within Octave itself for further |
| 508 details. | |
| 5164 | 509 |
|
8817
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
510 @node Return Types of Operators and Functions |
| 5164 | 511 @subsubsection The Return Types of Operators and Functions |
| 512 | |
| 5506 | 513 The two basic reasons to use sparse matrices are to reduce the memory |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
514 usage and to not have to do calculations on zero elements. The two are |
| 5164 | 515 closely related in that the computation time on a sparse matrix operator |
| 5506 | 516 or function is roughly linear with the number of non-zero elements. |
| 5164 | 517 |
| 518 Therefore, there is a certain density of non-zero elements of a matrix | |
| 519 where it no longer makes sense to store it as a sparse matrix, but rather | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
520 as a full matrix. For this reason operators and functions that have a |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
521 high probability of returning a full matrix will always return one. For |
| 5164 | 522 example adding a scalar constant to a sparse matrix will almost always |
| 523 make it a full matrix, and so the example | |
| 524 | |
| 525 @example | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
526 @group |
| 5164 | 527 speye(3) + 0 |
| 528 @result{} 1 0 0 | |
| 529 0 1 0 | |
| 530 0 0 1 | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
531 @end group |
| 5164 | 532 @end example |
| 533 | |
| 7330 | 534 returns a full matrix as can be seen. |
| 535 | |
| 536 | |
| 537 Additionally, if @code{sparse_auto_mutate} is true, all sparse functions | |
| 538 test the amount of memory occupied by the sparse matrix to see if the | |
| 539 amount of storage used is larger than the amount used by the full | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
540 equivalent. Therefore @code{speye (2) * 1} will return a full matrix as |
| 5164 | 541 the memory used is smaller for the full version than the sparse version. |
| 542 | |
| 543 As all of the mixed operators and functions between full and sparse | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
544 matrices exist, in general this does not cause any problems. However, |
| 5164 | 545 one area where it does cause a problem is where a sparse matrix is |
| 546 promoted to a full matrix, where subsequent operations would resparsify | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
547 the matrix. Such cases are rare, but can be artificially created, for |
| 5164 | 548 example @code{(fliplr(speye(3)) + speye(3)) - speye(3)} gives a full |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
549 matrix when it should give a sparse one. In general, where such cases |
| 5164 | 550 occur, they impose only a small memory penalty. |
| 551 | |
| 5648 | 552 There is however one known case where this behavior of Octave's |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
553 sparse matrices will cause a problem. That is in the handling of the |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
554 @dfn{diag} function. Whether @dfn{diag} returns a sparse or full matrix |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
555 depending on the type of its input arguments. So |
| 5164 | 556 |
| 557 @example | |
| 558 a = diag (sparse([1,2,3]), -1); | |
| 559 @end example | |
| 560 | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
561 should return a sparse matrix. To ensure this actually happens, the |
| 5164 | 562 @dfn{sparse} function, and other functions based on it like @dfn{speye}, |
| 563 always returns a sparse matrix, even if the memory used will be larger | |
| 564 than its full representation. | |
| 565 | |
| 7330 | 566 @DOCSTRING(sparse_auto_mutate) |
| 567 | |
| 568 Note that the @code{sparse_auto_mutate} option is incompatible with | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
569 @sc{matlab}, and so it is off by default. |
| 7330 | 570 |
|
8817
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
571 @node Mathematical Considerations |
| 5164 | 572 @subsubsection Mathematical Considerations |
| 573 | |
| 574 The attempt has been made to make sparse matrices behave in exactly the | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
575 same manner as there full counterparts. However, there are certain differences |
| 5164 | 576 and especially differences with other products sparse implementations. |
| 577 | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
578 Firstly, the "./" and ".^" operators must be used with care. Consider what |
| 5164 | 579 the examples |
| 580 | |
| 581 @example | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
582 @group |
| 5164 | 583 s = speye (4); |
| 584 a1 = s .^ 2; | |
| 585 a2 = s .^ s; | |
| 586 a3 = s .^ -2; | |
| 587 a4 = s ./ 2; | |
| 588 a5 = 2 ./ s; | |
| 589 a6 = s ./ s; | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
590 @end group |
| 5164 | 591 @end example |
| 592 | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
593 will give. The first example of @var{s} raised to the power of 2 causes |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
594 no problems. However @var{s} raised element-wise to itself involves a |
| 6431 | 595 large number of terms @code{0 .^ 0} which is 1. There @code{@var{s} .^ |
| 5164 | 596 @var{s}} is a full matrix. |
| 597 | |
|
8347
fa78cb8d8a5c
corrections for typos
Brian Gough<bjg@network-theory.co.uk>
parents:
8286
diff
changeset
|
598 Likewise @code{@var{s} .^ -2} involves terms like @code{0 .^ -2} which |
| 5164 | 599 is infinity, and so @code{@var{s} .^ -2} is equally a full matrix. |
| 600 | |
| 601 For the "./" operator @code{@var{s} ./ 2} has no problems, but | |
| 602 @code{2 ./ @var{s}} involves a large number of infinity terms as well | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
603 and is equally a full matrix. The case of @code{@var{s} ./ @var{s}} |
| 5164 | 604 involves terms like @code{0 ./ 0} which is a @code{NaN} and so this |
| 605 is equally a full matrix with the zero elements of @var{s} filled with | |
| 606 @code{NaN} values. | |
| 607 | |
| 5648 | 608 The above behavior is consistent with full matrices, but is not |
| 5164 | 609 consistent with sparse implementations in other products. |
| 610 | |
| 611 A particular problem of sparse matrices comes about due to the fact that | |
| 612 as the zeros are not stored, the sign-bit of these zeros is equally not | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
613 stored. In certain cases the sign-bit of zero is important. For example |
| 5164 | 614 |
| 615 @example | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
616 @group |
| 5164 | 617 a = 0 ./ [-1, 1; 1, -1]; |
| 618 b = 1 ./ a | |
| 619 @result{} -Inf Inf | |
| 620 Inf -Inf | |
| 621 c = 1 ./ sparse (a) | |
| 622 @result{} Inf Inf | |
| 623 Inf Inf | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
624 @end group |
| 5164 | 625 @end example |
| 626 | |
| 5648 | 627 To correct this behavior would mean that zero elements with a negative |
| 5164 | 628 sign-bit would need to be stored in the matrix to ensure that their |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
629 sign-bit was respected. This is not done at this time, for reasons of |
| 6750 | 630 efficiency, and so the user is warned that calculations where the sign-bit |
| 5164 | 631 of zero is important must not be done using sparse matrices. |
| 632 | |
| 5648 | 633 In general any function or operator used on a sparse matrix will |
| 634 result in a sparse matrix with the same or a larger number of non-zero | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
635 elements than the original matrix. This is particularly true for the |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
636 important case of sparse matrix factorizations. The usual way to |
| 5648 | 637 address this is to reorder the matrix, such that its factorization is |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
638 sparser than the factorization of the original matrix. That is the |
| 5648 | 639 factorization of @code{L * U = P * S * Q} has sparser terms @code{L} |
| 640 and @code{U} than the equivalent factorization @code{L * U = S}. | |
| 641 | |
| 642 Several functions are available to reorder depending on the type of the | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
643 matrix to be factorized. If the matrix is symmetric positive-definite, |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
644 then @dfn{symamd} or @dfn{csymamd} should be used. Otherwise |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
645 @dfn{amd}, @dfn{colamd} or @dfn{ccolamd} should be used. For completeness |
| 5648 | 646 the reordering functions @dfn{colperm} and @dfn{randperm} are |
| 647 also available. | |
| 648 | |
| 8829 | 649 @xref{fig:simplematrix}, for an example of the structure of a simple |
| 5648 | 650 positive definite matrix. |
| 5506 | 651 |
| 5648 | 652 @float Figure,fig:simplematrix |
|
9088
77e71f3da3d6
Fix documentation image printing under new development code
Rik <rdrider0-list@yahoo.com>
parents:
9066
diff
changeset
|
653 @center @image{spmatrix,4in} |
| 5648 | 654 @caption{Structure of simple sparse matrix.} |
| 655 @end float | |
| 5506 | 656 |
| 7001 | 657 The standard Cholesky factorization of this matrix can be |
| 5648 | 658 obtained by the same command that would be used for a full |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
659 matrix. This can be visualized with the command |
| 5652 | 660 @code{r = chol(A); spy(r);}. |
| 661 @ifset HAVE_CHOLMOD | |
| 662 @ifset HAVE_COLAMD | |
| 663 @xref{fig:simplechol}. | |
| 664 @end ifset | |
| 665 @end ifset | |
| 666 The original matrix had | |
| 5648 | 667 @ifinfo |
| 668 @ifnothtml | |
| 669 43 | |
| 670 @end ifnothtml | |
| 671 @end ifinfo | |
| 672 @ifset htmltex | |
| 673 598 | |
| 674 @end ifset | |
| 675 non-zero terms, while this Cholesky factorization has | |
| 676 @ifinfo | |
| 677 @ifnothtml | |
| 678 71, | |
| 679 @end ifnothtml | |
| 680 @end ifinfo | |
| 681 @ifset htmltex | |
| 682 10200, | |
| 683 @end ifset | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
684 with only half of the symmetric matrix being stored. This |
| 5648 | 685 is a significant level of fill in, and although not an issue |
| 686 for such a small test case, can represents a large overhead | |
| 687 in working with other sparse matrices. | |
| 5164 | 688 |
| 5648 | 689 The appropriate sparsity preserving permutation of the original |
| 690 matrix is given by @dfn{symamd} and the factorization using this | |
| 691 reordering can be visualized using the command @code{q = symamd(A); | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
692 r = chol(A(q,q)); spy(r)}. This gives |
| 5648 | 693 @ifinfo |
| 694 @ifnothtml | |
| 695 29 | |
| 696 @end ifnothtml | |
| 697 @end ifinfo | |
| 698 @ifset htmltex | |
| 699 399 | |
| 700 @end ifset | |
| 701 non-zero terms which is a significant improvement. | |
| 5164 | 702 |
| 5648 | 703 The Cholesky factorization itself can be used to determine the |
| 704 appropriate sparsity preserving reordering of the matrix during the | |
| 705 factorization, In that case this might be obtained with three return | |
| 706 arguments as r@code{[r, p, q] = chol(A); spy(r)}. | |
| 5164 | 707 |
| 5648 | 708 @ifset HAVE_CHOLMOD |
| 709 @ifset HAVE_COLAMD | |
| 710 @float Figure,fig:simplechol | |
|
9088
77e71f3da3d6
Fix documentation image printing under new development code
Rik <rdrider0-list@yahoo.com>
parents:
9066
diff
changeset
|
711 @center @image{spchol,4in} |
| 5648 | 712 @caption{Structure of the un-permuted Cholesky factorization of the above matrix.} |
| 713 @end float | |
| 5164 | 714 |
| 5648 | 715 @float Figure,fig:simplecholperm |
|
9088
77e71f3da3d6
Fix documentation image printing under new development code
Rik <rdrider0-list@yahoo.com>
parents:
9066
diff
changeset
|
716 @center @image{spcholperm,4in} |
| 5648 | 717 @caption{Structure of the permuted Cholesky factorization of the above matrix.} |
| 718 @end float | |
| 719 @end ifset | |
| 720 @end ifset | |
| 5164 | 721 |
| 5648 | 722 In the case of an asymmetric matrix, the appropriate sparsity |
| 723 preserving permutation is @dfn{colamd} and the factorization using | |
| 724 this reordering can be visualized using the command @code{q = | |
| 725 colamd(A); [l, u, p] = lu(A(:,q)); spy(l+u)}. | |
| 5164 | 726 |
| 5648 | 727 Finally, Octave implicitly reorders the matrix when using the div (/) |
| 728 and ldiv (\) operators, and so no the user does not need to explicitly | |
| 729 reorder the matrix to maximize performance. | |
| 730 | |
| 7619 | 731 @DOCSTRING(amd) |
| 732 | |
| 6620 | 733 @DOCSTRING(ccolamd) |
| 734 | |
| 735 @DOCSTRING(colamd) | |
| 736 | |
| 737 @DOCSTRING(colperm) | |
| 738 | |
| 739 @DOCSTRING(csymamd) | |
| 740 | |
| 741 @DOCSTRING(dmperm) | |
| 742 | |
| 743 @DOCSTRING(symamd) | |
| 744 | |
| 745 @DOCSTRING(symrcm) | |
| 746 | |
|
8817
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
747 @node Sparse Linear Algebra |
| 5164 | 748 @section Linear Algebra on Sparse Matrices |
| 749 | |
|
8488
cdb4788879b3
[docs] poly-morphic => polymorphic
Brian Gough <bjg@gnu.org>
parents:
8417
diff
changeset
|
750 Octave includes a polymorphic solver for sparse matrices, where |
| 5164 | 751 the exact solver used to factorize the matrix, depends on the properties |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
752 of the sparse matrix itself. Generally, the cost of determining the matrix type |
| 5322 | 753 is small relative to the cost of factorizing the matrix itself, but in any |
| 754 case the matrix type is cached once it is calculated, so that it is not | |
| 755 re-determined each time it is used in a linear equation. | |
| 5164 | 756 |
| 757 The selection tree for how the linear equation is solve is | |
| 758 | |
| 759 @enumerate 1 | |
| 5648 | 760 @item If the matrix is diagonal, solve directly and goto 8 |
| 5164 | 761 |
| 762 @item If the matrix is a permuted diagonal, solve directly taking into | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
763 account the permutations. Goto 8 |
| 5164 | 764 |
| 5648 | 765 @item If the matrix is square, banded and if the band density is less |
| 766 than that given by @code{spparms ("bandden")} continue, else goto 4. | |
| 5164 | 767 |
| 768 @enumerate a | |
| 769 @item If the matrix is tridiagonal and the right-hand side is not sparse | |
| 5648 | 770 continue, else goto 3b. |
| 5164 | 771 |
| 772 @enumerate | |
| 773 @item If the matrix is hermitian, with a positive real diagonal, attempt | |
|
9088
77e71f3da3d6
Fix documentation image printing under new development code
Rik <rdrider0-list@yahoo.com>
parents:
9066
diff
changeset
|
774 Cholesky factorization using @sc{lapack} xPTSV. |
| 5164 | 775 |
| 776 @item If the above failed or the matrix is not hermitian with a positive | |
| 777 real diagonal use Gaussian elimination with pivoting using | |
|
9088
77e71f3da3d6
Fix documentation image printing under new development code
Rik <rdrider0-list@yahoo.com>
parents:
9066
diff
changeset
|
778 @sc{lapack} xGTSV, and goto 8. |
| 5164 | 779 @end enumerate |
| 780 | |
| 781 @item If the matrix is hermitian with a positive real diagonal, attempt | |
|
9088
77e71f3da3d6
Fix documentation image printing under new development code
Rik <rdrider0-list@yahoo.com>
parents:
9066
diff
changeset
|
782 Cholesky factorization using @sc{lapack} xPBTRF. |
| 5164 | 783 |
| 784 @item if the above failed or the matrix is not hermitian with a positive | |
| 785 real diagonal use Gaussian elimination with pivoting using | |
|
9088
77e71f3da3d6
Fix documentation image printing under new development code
Rik <rdrider0-list@yahoo.com>
parents:
9066
diff
changeset
|
786 @sc{lapack} xGBTRF, and goto 8. |
| 5164 | 787 @end enumerate |
| 788 | |
| 789 @item If the matrix is upper or lower triangular perform a sparse forward | |
| 5648 | 790 or backward substitution, and goto 8 |
| 5164 | 791 |
| 5322 | 792 @item If the matrix is a upper triangular matrix with column permutations |
| 793 or lower triangular matrix with row permutations, perform a sparse forward | |
| 5648 | 794 or backward substitution, and goto 8 |
| 5164 | 795 |
| 5648 | 796 @item If the matrix is square, hermitian with a real positive diagonal, attempt |
| 5506 | 797 sparse Cholesky factorization using CHOLMOD. |
| 5164 | 798 |
| 799 @item If the sparse Cholesky factorization failed or the matrix is not | |
| 5648 | 800 hermitian with a real positive diagonal, and the matrix is square, factorize |
| 801 using UMFPACK. | |
| 5164 | 802 |
| 803 @item If the matrix is not square, or any of the previous solvers flags | |
| 5648 | 804 a singular or near singular matrix, find a minimum norm solution using |
| 7096 | 805 CXSPARSE@footnote{The CHOLMOD, UMFPACK and CXSPARSE packages were |
| 806 written by Tim Davis and are available at | |
| 807 http://www.cise.ufl.edu/research/sparse/}. | |
| 5164 | 808 @end enumerate |
| 809 | |
| 810 The band density is defined as the number of non-zero values in the matrix | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
811 divided by the number of non-zero values in the matrix. The banded matrix |
| 5164 | 812 solvers can be entirely disabled by using @dfn{spparms} to set @code{bandden} |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
813 to 1 (i.e., @code{spparms ("bandden", 1)}). |
| 5164 | 814 |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
815 The QR solver factorizes the problem with a Dulmage-Mendelsohn, to |
| 6939 | 816 separate the problem into blocks that can be treated as over-determined, |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
817 multiple well determined blocks, and a final over-determined block. For |
| 6939 | 818 matrices with blocks of strongly connected nodes this is a big win as |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
819 LU decomposition can be used for many blocks. It also significantly |
| 5681 | 820 improves the chance of finding a solution to over-determined problems |
| 821 rather than just returning a vector of @dfn{NaN}'s. | |
| 822 | |
| 823 All of the solvers above, can calculate an estimate of the condition | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
824 number. This can be used to detect numerical stability problems in the |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
825 solution and force a minimum norm solution to be used. However, for |
| 5681 | 826 narrow banded, triangular or diagonal matrices, the cost of |
| 827 calculating the condition number is significant, and can in fact | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
828 exceed the cost of factoring the matrix. Therefore the condition |
| 6939 | 829 number is not calculated in these cases, and Octave relies on simpler |
|
9088
77e71f3da3d6
Fix documentation image printing under new development code
Rik <rdrider0-list@yahoo.com>
parents:
9066
diff
changeset
|
830 techniques to detect singular matrices or the underlying @sc{lapack} code in |
| 5681 | 831 the case of banded matrices. |
| 5164 | 832 |
| 5322 | 833 The user can force the type of the matrix with the @code{matrix_type} |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
834 function. This overcomes the cost of discovering the type of the matrix. |
|
8347
fa78cb8d8a5c
corrections for typos
Brian Gough<bjg@network-theory.co.uk>
parents:
8286
diff
changeset
|
835 However, it should be noted that identifying the type of the matrix incorrectly |
| 5322 | 836 will lead to unpredictable results, and so @code{matrix_type} should be |
| 5506 | 837 used with care. |
| 5322 | 838 |
| 6620 | 839 @DOCSTRING(normest) |
| 840 | |
|
8286
6f2d95255911
fix @seealso references to point to existing anchors
Thorsten Meyer <thorsten.meyier@gmx.de>
parents:
8106
diff
changeset
|
841 @DOCSTRING(onenormest) |
|
6f2d95255911
fix @seealso references to point to existing anchors
Thorsten Meyer <thorsten.meyier@gmx.de>
parents:
8106
diff
changeset
|
842 |
| 7189 | 843 @DOCSTRING(condest) |
| 844 | |
| 6620 | 845 @DOCSTRING(spparms) |
| 846 | |
| 847 @DOCSTRING(sprank) | |
| 848 | |
| 849 @DOCSTRING(symbfact) | |
| 850 | |
|
7681
b1c1133641ee
Add the spaugment function
David Bateman <dbateman@free.fr>
parents:
7619
diff
changeset
|
851 For non square matrices, the user can also utilize the @code{spaugment} |
|
b1c1133641ee
Add the spaugment function
David Bateman <dbateman@free.fr>
parents:
7619
diff
changeset
|
852 function to find a least squares solution to a linear equation. |
|
b1c1133641ee
Add the spaugment function
David Bateman <dbateman@free.fr>
parents:
7619
diff
changeset
|
853 |
|
b1c1133641ee
Add the spaugment function
David Bateman <dbateman@free.fr>
parents:
7619
diff
changeset
|
854 @DOCSTRING(spaugment) |
|
b1c1133641ee
Add the spaugment function
David Bateman <dbateman@free.fr>
parents:
7619
diff
changeset
|
855 |
|
8417
654bcfb937bf
Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
8347
diff
changeset
|
856 Finally, the function @code{eigs} can be used to calculate a limited |
|
654bcfb937bf
Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
8347
diff
changeset
|
857 number of eigenvalues and eigenvectors based on a selection criteria |
|
654bcfb937bf
Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
8347
diff
changeset
|
858 and likewise for @code{svds} which calculates a limited number of |
|
654bcfb937bf
Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
8347
diff
changeset
|
859 singular values and vectors. |
|
654bcfb937bf
Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
8347
diff
changeset
|
860 |
|
654bcfb937bf
Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
8347
diff
changeset
|
861 @DOCSTRING(eigs) |
|
654bcfb937bf
Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
8347
diff
changeset
|
862 |
|
654bcfb937bf
Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
8347
diff
changeset
|
863 @DOCSTRING(svds) |
|
654bcfb937bf
Add the eigs and svds functions
David Bateman <dbateman@free.fr>
parents:
8347
diff
changeset
|
864 |
|
8817
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
865 @node Iterative Techniques |
| 5164 | 866 @section Iterative Techniques applied to sparse matrices |
| 867 | |
| 6620 | 868 The left division @code{\} and right division @code{/} operators, |
| 869 discussed in the previous section, use direct solvers to resolve a | |
| 870 linear equation of the form @code{@var{x} = @var{A} \ @var{b}} or | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
871 @code{@var{x} = @var{b} / @var{A}}. Octave equally includes a number of |
| 6620 | 872 functions to solve sparse linear equations using iterative techniques. |
| 873 | |
| 874 @DOCSTRING(pcg) | |
| 875 | |
| 876 @DOCSTRING(pcr) | |
| 5837 | 877 |
| 6620 | 878 The speed with which an iterative solver converges to a solution can be |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
879 accelerated with the use of a pre-conditioning matrix @var{M}. In this |
| 6620 | 880 case the linear equation @code{@var{M}^-1 * @var{x} = @var{M}^-1 * |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
881 @var{A} \ @var{b}} is solved instead. Typical pre-conditioning matrices |
| 6620 | 882 are partial factorizations of the original matrix. |
| 5648 | 883 |
| 6620 | 884 @DOCSTRING(luinc) |
| 885 | |
|
8817
03b7f618ab3d
include docstrings for new functions in the manual
John W. Eaton <jwe@octave.org>
parents:
8488
diff
changeset
|
886 @node Real Life Example |
| 5648 | 887 @section Real Life Example of the use of Sparse Matrices |
| 888 | |
| 889 A common application for sparse matrices is in the solution of Finite | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
890 Element Models. Finite element models allow numerical solution of |
| 5648 | 891 partial differential equations that do not have closed form solutions, |
| 892 typically because of the complex shape of the domain. | |
| 893 | |
| 894 In order to motivate this application, we consider the boundary value | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
895 Laplace equation. This system can model scalar potential fields, such |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
896 as heat or electrical potential. Given a medium |
| 5648 | 897 @iftex |
| 898 @tex | |
| 899 $\Omega$ | |
| 900 @end tex | |
| 901 @end iftex | |
| 902 @ifinfo | |
| 903 Omega | |
| 904 @end ifinfo | |
| 905 with boundary | |
| 906 @iftex | |
| 907 @tex | |
| 908 $\partial\Omega$ | |
| 909 @end tex | |
| 910 @end iftex | |
| 911 @ifinfo | |
| 912 dOmega | |
| 913 @end ifinfo | |
| 914 . At all points on the | |
| 915 @iftex | |
| 916 @tex | |
| 917 $\partial\Omega$ | |
| 918 @end tex | |
| 919 @end iftex | |
| 920 @ifinfo | |
| 921 dOmega | |
| 922 @end ifinfo | |
| 923 the boundary conditions are known, and we wish to calculate the potential in | |
| 924 @iftex | |
| 925 @tex | |
| 926 $\Omega$ | |
| 927 @end tex | |
| 928 @end iftex | |
| 929 @ifinfo | |
| 930 Omega | |
| 931 @end ifinfo | |
| 932 . Boundary conditions may specify the potential (Dirichlet | |
| 933 boundary condition), its normal derivative across the boundary | |
| 934 (Neumann boundary condition), or a weighted sum of the potential and | |
| 935 its derivative (Cauchy boundary condition). | |
| 936 | |
| 937 In a thermal model, we want to calculate the temperature in | |
| 938 @iftex | |
| 939 @tex | |
| 940 $\Omega$ | |
| 941 @end tex | |
| 942 @end iftex | |
| 943 @ifinfo | |
| 944 Omega | |
| 945 @end ifinfo | |
| 946 and know the boundary temperature (Dirichlet condition) | |
| 947 or heat flux (from which we can calculate the Neumann condition | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
948 by dividing by the thermal conductivity at the boundary). Similarly, |
| 5648 | 949 in an electrical model, we want to calculate the voltage in |
| 950 @iftex | |
| 951 @tex | |
| 952 $\Omega$ | |
| 953 @end tex | |
| 954 @end iftex | |
| 955 @ifinfo | |
| 956 Omega | |
| 957 @end ifinfo | |
| 958 and know the boundary voltage (Dirichlet) or current | |
| 959 (Neumann condition after diving by the electrical conductivity). | |
| 960 In an electrical model, it is common for much of the boundary | |
| 961 to be electrically isolated; this is a Neumann boundary condition | |
| 962 with the current equal to zero. | |
| 963 | |
| 964 The simplest finite element models will divide | |
| 965 @iftex | |
| 966 @tex | |
| 967 $\Omega$ | |
| 968 @end tex | |
| 969 @end iftex | |
| 970 @ifinfo | |
| 971 Omega | |
| 972 @end ifinfo | |
| 973 into simplexes (triangles in 2D, pyramids in 3D). | |
| 974 @ifset htmltex | |
| 975 We take as an 3D example a cylindrical liquid filled tank with a small | |
| 976 non-conductive ball from the EIDORS project@footnote{EIDORS - Electrical | |
| 977 Impedance Tomography and Diffuse optical Tomography Reconstruction Software | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
978 @url{http://eidors3d.sourceforge.net}}. This is model is designed to reflect |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
979 an application of electrical impedance tomography, where current patterns |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
980 are applied to such a tank in order to image the internal conductivity |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
981 distribution. In order to describe the FEM geometry, we have a matrix of |
| 5648 | 982 vertices @code{nodes} and simplices @code{elems}. |
| 983 @end ifset | |
| 984 | |
| 985 The following example creates a simple rectangular 2D electrically | |
| 986 conductive medium with 10 V and 20 V imposed on opposite sides | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
987 (Dirichlet boundary conditions). All other edges are electrically |
| 5648 | 988 isolated. |
| 989 | |
| 990 @example | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
991 @group |
| 5648 | 992 node_y= [1;1.2;1.5;1.8;2]*ones(1,11); |
| 993 node_x= ones(5,1)*[1,1.05,1.1,1.2, ... | |
| 994 1.3,1.5,1.7,1.8,1.9,1.95,2]; | |
| 995 nodes= [node_x(:), node_y(:)]; | |
| 996 | |
| 997 [h,w]= size(node_x); | |
| 998 elems= []; | |
| 999 for idx= 1:w-1 | |
| 1000 widx= (idx-1)*h; | |
| 1001 elems= [elems; ... | |
| 1002 widx+[(1:h-1);(2:h);h+(1:h-1)]'; ... | |
| 1003 widx+[(2:h);h+(2:h);h+(1:h-1)]' ]; | |
| 1004 endfor | |
| 1005 | |
| 1006 E= size(elems,1); # No. of simplices | |
| 1007 N= size(nodes,1); # No. of vertices | |
| 1008 D= size(elems,2); # dimensions+1 | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1009 @end group |
| 5648 | 1010 @end example |
| 1011 | |
| 1012 This creates a N-by-2 matrix @code{nodes} and a E-by-3 matrix | |
| 1013 @code{elems} with values, which define finite element triangles: | |
| 5164 | 1014 |
| 5648 | 1015 @example |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1016 @group |
| 5648 | 1017 nodes(1:7,:)' |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1018 1.00 1.00 1.00 1.00 1.00 1.05 1.05 @dots{} |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1019 1.00 1.20 1.50 1.80 2.00 1.00 1.20 @dots{} |
| 5648 | 1020 |
| 1021 elems(1:7,:)' | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1022 1 2 3 4 2 3 4 @dots{} |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1023 2 3 4 5 7 8 9 @dots{} |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1024 6 7 8 9 6 7 8 @dots{} |
|
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1025 @end group |
| 5648 | 1026 @end example |
| 1027 | |
| 1028 Using a first order FEM, we approximate the electrical conductivity | |
| 1029 distribution in | |
| 1030 @iftex | |
| 1031 @tex | |
| 1032 $\Omega$ | |
| 1033 @end tex | |
| 1034 @end iftex | |
| 1035 @ifinfo | |
| 1036 Omega | |
| 1037 @end ifinfo | |
| 1038 as constant on each simplex (represented by the vector @code{conductivity}). | |
| 1039 Based on the finite element geometry, we first calculate a system (or | |
| 1040 stiffness) matrix for each simplex (represented as 3-by-3 elements on the | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1041 diagonal of the element-wise system matrix @code{SE}. Based on @code{SE} |
| 5648 | 1042 and a N-by-DE connectivity matrix @code{C}, representing the connections |
| 7001 | 1043 between simplices and vertices, the global connectivity matrix @code{S} is |
| 5648 | 1044 calculated. |
| 1045 | |
| 1046 @example | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1047 @group |
| 5648 | 1048 # Element conductivity |
| 1049 conductivity= [1*ones(1,16), ... | |
| 1050 2*ones(1,48), 1*ones(1,16)]; | |
| 1051 | |
| 1052 # Connectivity matrix | |
| 1053 C = sparse ((1:D*E), reshape (elems', ... | |
| 1054 D*E, 1), 1, D*E, N); | |
| 1055 | |
| 1056 # Calculate system matrix | |
| 1057 Siidx = floor ([0:D*E-1]'/D) * D * ... | |
| 1058 ones(1,D) + ones(D*E,1)*(1:D) ; | |
| 1059 Sjidx = [1:D*E]'*ones(1,D); | |
| 1060 Sdata = zeros(D*E,D); | |
| 1061 dfact = factorial(D-1); | |
| 1062 for j=1:E | |
| 1063 a = inv([ones(D,1), ... | |
| 1064 nodes(elems(j,:), :)]); | |
| 1065 const = conductivity(j) * 2 / ... | |
| 1066 dfact / abs(det(a)); | |
| 1067 Sdata(D*(j-1)+(1:D),:) = const * ... | |
| 1068 a(2:D,:)' * a(2:D,:); | |
| 1069 endfor | |
| 1070 # Element-wise system matrix | |
| 1071 SE= sparse(Siidx,Sjidx,Sdata); | |
| 1072 # Global system matrix | |
| 1073 S= C'* SE *C; | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1074 @end group |
| 5648 | 1075 @end example |
| 1076 | |
| 1077 The system matrix acts like the conductivity | |
| 1078 @iftex | |
| 1079 @tex | |
| 1080 $S$ | |
| 1081 @end tex | |
| 1082 @end iftex | |
| 1083 @ifinfo | |
| 1084 @code{S} | |
| 1085 @end ifinfo | |
| 1086 in Ohm's law | |
| 1087 @iftex | |
| 1088 @tex | |
| 1089 $SV = I$. | |
| 1090 @end tex | |
| 1091 @end iftex | |
| 1092 @ifinfo | |
| 1093 @code{S * V = I}. | |
| 1094 @end ifinfo | |
| 1095 Based on the Dirichlet and Neumann boundary conditions, we are able to | |
| 1096 solve for the voltages at each vertex @code{V}. | |
| 1097 | |
| 1098 @example | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1099 @group |
| 5648 | 1100 # Dirichlet boundary conditions |
| 1101 D_nodes=[1:5, 51:55]; | |
| 1102 D_value=[10*ones(1,5), 20*ones(1,5)]; | |
| 1103 | |
| 1104 V= zeros(N,1); | |
| 1105 V(D_nodes) = D_value; | |
| 1106 idx = 1:N; # vertices without Dirichlet | |
| 1107 # boundary condns | |
| 1108 idx(D_nodes) = []; | |
| 1109 | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1110 # Neumann boundary conditions. Note that |
| 5648 | 1111 # N_value must be normalized by the |
| 1112 # boundary length and element conductivity | |
| 1113 N_nodes=[]; | |
| 1114 N_value=[]; | |
| 1115 | |
| 1116 Q = zeros(N,1); | |
| 1117 Q(N_nodes) = N_value; | |
| 1118 | |
| 1119 V(idx) = S(idx,idx) \ ( Q(idx) - ... | |
| 1120 S(idx,D_nodes) * V(D_nodes)); | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1121 @end group |
| 5648 | 1122 @end example |
| 1123 | |
| 1124 Finally, in order to display the solution, we show each solved voltage | |
| 1125 value in the z-axis for each simplex vertex. | |
| 1126 @ifset htmltex | |
| 5652 | 1127 @ifset HAVE_CHOLMOD |
| 1128 @ifset HAVE_UMFPACK | |
| 1129 @ifset HAVE_COLAMD | |
| 5648 | 1130 @xref{fig:femmodel}. |
| 1131 @end ifset | |
| 5652 | 1132 @end ifset |
| 1133 @end ifset | |
| 1134 @end ifset | |
| 5648 | 1135 |
| 1136 @example | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1137 @group |
| 5648 | 1138 elemx = elems(:,[1,2,3,1])'; |
| 1139 xelems = reshape (nodes(elemx, 1), 4, E); | |
| 1140 yelems = reshape (nodes(elemx, 2), 4, E); | |
| 1141 velems = reshape (V(elemx), 4, E); | |
| 1142 plot3 (xelems,yelems,velems,'k'); | |
| 1143 print ('grid.eps'); | |
|
9066
be150a172010
Cleanup documentation for diagperm.texi, sparse.texi
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
1144 @end group |
| 5648 | 1145 @end example |
| 1146 | |
| 1147 | |
| 1148 @ifset htmltex | |
| 1149 @ifset HAVE_CHOLMOD | |
| 1150 @ifset HAVE_UMFPACK | |
| 1151 @ifset HAVE_COLAMD | |
| 1152 @float Figure,fig:femmodel | |
|
9088
77e71f3da3d6
Fix documentation image printing under new development code
Rik <rdrider0-list@yahoo.com>
parents:
9066
diff
changeset
|
1153 @center @image{grid,4in} |
| 5648 | 1154 @caption{Example finite element model the showing triangular elements. |
| 1155 The height of each vertex corresponds to the solution value.} | |
| 1156 @end float | |
| 1157 @end ifset | |
| 1158 @end ifset | |
| 1159 @end ifset | |
| 1160 @end ifset |