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