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