Mercurial > octave-nkf

comparison doc/interpreter/sparse.txi @ 19627:446c46af4b42 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
strip trailing whitespace from most source files * Makefile.am, NEWS, build-aux/common.mk, configure.ac, doc/Makefile.am, doc/doxyhtml/Makefile.am, doc/interpreter/Makefile.am, doc/interpreter/arith.txi, doc/interpreter/audio.txi, doc/interpreter/basics.txi, doc/interpreter/bugs.txi, doc/interpreter/container.txi, doc/interpreter/cp-idx.txi, doc/interpreter/data.txi, doc/interpreter/debug.txi, doc/interpreter/diagperm.txi, doc/interpreter/diffeq.txi, doc/interpreter/doccheck/README, doc/interpreter/doccheck/spellcheck, doc/interpreter/emacs.txi, doc/interpreter/errors.txi, doc/interpreter/eval.txi, doc/interpreter/expr.txi, doc/interpreter/external.txi, doc/interpreter/fn-idx.txi, doc/interpreter/func.txi, doc/interpreter/geometry.txi, doc/interpreter/geometryimages.m, doc/interpreter/gpl.txi, doc/interpreter/grammar.txi, doc/interpreter/gui.txi, doc/interpreter/image.txi, doc/interpreter/install.txi, doc/interpreter/interp.txi, doc/interpreter/interpimages.m, doc/interpreter/intro.txi, doc/interpreter/io.txi, doc/interpreter/java.txi, doc/interpreter/linalg.txi, doc/interpreter/macros.texi, doc/interpreter/matrix.txi, doc/interpreter/munge-texi.pl, doc/interpreter/nonlin.txi, doc/interpreter/numbers.txi, doc/interpreter/obsolete.txi, doc/interpreter/octave-config.1, doc/interpreter/octave.texi, doc/interpreter/oop.txi, doc/interpreter/op-idx.txi, doc/interpreter/optim.txi, doc/interpreter/package.txi, doc/interpreter/plot.txi, doc/interpreter/poly.txi, doc/interpreter/preface.txi, doc/interpreter/quad.txi, doc/interpreter/set.txi, doc/interpreter/signal.txi, doc/interpreter/sparse.txi, doc/interpreter/sparseimages.m, doc/interpreter/splineimages.m, doc/interpreter/stats.txi, doc/interpreter/stmt.txi, doc/interpreter/strings.txi, doc/interpreter/system.txi, doc/interpreter/testfun.txi, doc/interpreter/tips.txi, doc/interpreter/var.txi, doc/interpreter/vectorize.txi, doc/liboctave/Makefile.am, doc/liboctave/array.texi, doc/liboctave/bugs.texi, doc/liboctave/cp-idx.texi, doc/liboctave/dae.texi, doc/liboctave/diffeq.texi, doc/liboctave/error.texi, doc/liboctave/factor.texi, doc/liboctave/fn-idx.texi, doc/liboctave/gpl.texi, doc/liboctave/install.texi, doc/liboctave/intro.texi, doc/liboctave/liboctave.texi, doc/liboctave/matvec.texi, doc/liboctave/nleqn.texi, doc/liboctave/nlfunc.texi, doc/liboctave/ode.texi, doc/liboctave/optim.texi, doc/liboctave/preface.texi, doc/liboctave/quad.texi, doc/liboctave/range.texi, doc/refcard/Makefile.am, doc/refcard/refcard.tex, etc/HACKING, etc/NEWS.1, etc/NEWS.2, etc/NEWS.3, etc/OLD-ChangeLogs/ChangeLog, etc/OLD-ChangeLogs/doc-ChangeLog, etc/OLD-ChangeLogs/scripts-ChangeLog, etc/OLD-ChangeLogs/src-ChangeLog, etc/OLD-ChangeLogs/test-ChangeLog, etc/PROJECTS, etc/README.Cygwin, etc/README.MacOS, etc/README.MinGW, etc/README.gnuplot, etc/gdbinit, etc/icons/Makefile.am, examples/@polynomial/end.m, examples/@polynomial/subsasgn.m, examples/Makefile.am, examples/standalonebuiltin.cc, libgui/Makefile.am, libgui/qterminal/libqterminal/README, libgui/qterminal/libqterminal/unix/BlockArray.cpp, libgui/qterminal/libqterminal/unix/BlockArray.h, libgui/qterminal/libqterminal/unix/Character.h, libgui/qterminal/libqterminal/unix/CharacterColor.h, libgui/qterminal/libqterminal/unix/Emulation.cpp, libgui/qterminal/libqterminal/unix/Emulation.h, libgui/qterminal/libqterminal/unix/Filter.cpp, libgui/qterminal/libqterminal/unix/Filter.h, libgui/qterminal/libqterminal/unix/History.cpp, libgui/qterminal/libqterminal/unix/History.h, libgui/qterminal/libqterminal/unix/KeyboardTranslator.cpp, libgui/qterminal/libqterminal/unix/KeyboardTranslator.h, libgui/qterminal/libqterminal/unix/LineFont.h, libgui/qterminal/libqterminal/unix/QUnixTerminalImpl.cpp, libgui/qterminal/libqterminal/unix/QUnixTerminalImpl.h, libgui/qterminal/libqterminal/unix/Screen.cpp, libgui/qterminal/libqterminal/unix/Screen.h, libgui/qterminal/libqterminal/unix/ScreenWindow.cpp, libgui/qterminal/libqterminal/unix/ScreenWindow.h, libgui/qterminal/libqterminal/unix/TerminalCharacterDecoder.cpp, libgui/qterminal/libqterminal/unix/TerminalCharacterDecoder.h, libgui/qterminal/libqterminal/unix/Vt102Emulation.h, libgui/qterminal/libqterminal/win32/QWinTerminalImpl.cpp, libgui/qterminal/qterminal/main.cpp, libgui/src/m-editor/file-editor-tab.cc, libgui/src/octave-gui.cc, libgui/src/octave-qt-link.cc, libinterp/corefcn/data.cc, libinterp/corefcn/defun-int.h, libinterp/corefcn/det.cc, libinterp/corefcn/gl2ps-renderer.cc, libinterp/corefcn/graphics.cc, libinterp/corefcn/graphics.in.h, libinterp/corefcn/ls-mat5.cc, libinterp/corefcn/lu.cc, libinterp/corefcn/oct-tex-parser.yy, libinterp/corefcn/oct-tex-symbols.in, libinterp/corefcn/quadcc.cc, libinterp/corefcn/zfstream.cc, libinterp/dldfcn/__eigs__.cc, libinterp/dldfcn/__voronoi__.cc, libinterp/gendoc.pl, libinterp/genprops.awk, libinterp/mk-errno-list, libinterp/mk-pkg-add, libinterp/mkbuiltins, libinterp/mkdefs, libinterp/mkdocs, libinterp/mkops, libinterp/octave-value/ov-java.cc, libinterp/parse-tree/lex.ll, libinterp/parse-tree/oct-parse.in.yy, libinterp/parse-tree/octave.gperf, liboctave/Makefile.am, liboctave/array/Array.cc, liboctave/array/module.mk, liboctave/cruft/daspk/datv.f, liboctave/cruft/daspk/dcnst0.f, liboctave/cruft/daspk/dcnstr.f, liboctave/cruft/daspk/ddasic.f, liboctave/cruft/daspk/ddasid.f, liboctave/cruft/daspk/ddasik.f, liboctave/cruft/daspk/ddaspk.f, liboctave/cruft/daspk/ddstp.f, liboctave/cruft/daspk/ddwnrm.f, liboctave/cruft/daspk/dfnrmd.f, liboctave/cruft/daspk/dfnrmk.f, liboctave/cruft/daspk/dhels.f, liboctave/cruft/daspk/dheqr.f, liboctave/cruft/daspk/dinvwt.f, liboctave/cruft/daspk/dlinsd.f, liboctave/cruft/daspk/dlinsk.f, liboctave/cruft/daspk/dmatd.f, liboctave/cruft/daspk/dnedd.f, liboctave/cruft/daspk/dnedk.f, liboctave/cruft/daspk/dnsd.f, liboctave/cruft/daspk/dnsid.f, liboctave/cruft/daspk/dnsik.f, liboctave/cruft/daspk/dnsk.f, liboctave/cruft/daspk/dorth.f, liboctave/cruft/daspk/dslvd.f, liboctave/cruft/daspk/dslvk.f, liboctave/cruft/daspk/dspigm.f, liboctave/cruft/daspk/dyypnw.f, liboctave/cruft/dasrt/ddasrt.f, liboctave/cruft/dasrt/drchek.f, liboctave/cruft/dassl/ddaslv.f, liboctave/cruft/dassl/ddassl.f, liboctave/cruft/misc/blaswrap.c, liboctave/cruft/misc/module.mk, liboctave/cruft/odepack/cfode.f, liboctave/cruft/odepack/dlsode.f, liboctave/cruft/odepack/ewset.f, liboctave/cruft/odepack/intdy.f, liboctave/cruft/odepack/prepj.f, liboctave/cruft/odepack/sintdy.f, liboctave/cruft/odepack/slsode.f, liboctave/cruft/odepack/solsy.f, liboctave/cruft/odepack/ssolsy.f, liboctave/cruft/odepack/stode.f, liboctave/cruft/odepack/vnorm.f, liboctave/cruft/ranlib/Basegen.doc, liboctave/cruft/ranlib/README, liboctave/cruft/ranlib/genbet.f, liboctave/cruft/ranlib/genexp.f, liboctave/cruft/ranlib/gennch.f, liboctave/cruft/ranlib/gennf.f, liboctave/cruft/ranlib/gennor.f, liboctave/cruft/ranlib/getsd.f, liboctave/cruft/ranlib/initgn.f, liboctave/cruft/ranlib/phrtsd.f, liboctave/cruft/ranlib/randlib.fdoc, liboctave/cruft/ranlib/setsd.f, liboctave/cruft/ranlib/tstgmn.for, liboctave/cruft/ranlib/tstmid.for, liboctave/cruft/slatec-fn/atanh.f, liboctave/cruft/slatec-fn/datanh.f, liboctave/cruft/slatec-fn/xgmainc.f, liboctave/cruft/slatec-fn/xsgmainc.f, liboctave/numeric/module.mk, liboctave/operators/mk-ops.awk, liboctave/operators/mx-ops, liboctave/operators/sparse-mk-ops.awk, liboctave/operators/sparse-mx-ops, liboctave/operators/vx-ops, liboctave/util/module.mk, run-octave.in, scripts/@ftp/ftp.m, scripts/audio/wavread.m, scripts/deprecated/java_convert_matrix.m, scripts/deprecated/java_debug.m, scripts/deprecated/java_invoke.m, scripts/deprecated/java_new.m, scripts/deprecated/java_unsigned_conversion.m, scripts/deprecated/javafields.m, scripts/deprecated/javamethods.m, scripts/deprecated/shell_cmd.m, scripts/general/accumarray.m, scripts/general/display.m, scripts/general/fieldnames.m, scripts/general/interp1.m, scripts/general/interp2.m, scripts/general/interp3.m, scripts/general/isa.m, scripts/general/methods.m, scripts/general/sortrows.m, scripts/geometry/convhull.m, scripts/geometry/delaunay.m, scripts/geometry/delaunay3.m, scripts/geometry/delaunayn.m, scripts/geometry/griddata.m, scripts/geometry/griddatan.m, scripts/geometry/voronoi.m, scripts/geometry/voronoin.m, scripts/gui/guihandles.m, scripts/gui/inputdlg.m, scripts/gui/listdlg.m, scripts/gui/msgbox.m, scripts/gui/questdlg.m, scripts/gui/uigetfile.m, scripts/gui/waitbar.m, scripts/gui/warndlg.m, scripts/help/doc.m, scripts/help/help.m, scripts/help/type.m, scripts/image/bone.m, scripts/image/cmpermute.m, scripts/image/cmunique.m, scripts/image/colorcube.m, scripts/image/colormap.m, scripts/image/contrast.m, scripts/image/gray2ind.m, scripts/image/image.m, scripts/image/imshow.m, scripts/image/ind2gray.m, scripts/image/jet.m, scripts/image/rgb2ntsc.m, scripts/image/spinmap.m, scripts/io/importdata.m, scripts/io/strread.m, scripts/io/textread.m, scripts/io/textscan.m, scripts/java/java_get.m, scripts/java/java_set.m, scripts/java/javaaddpath.m, scripts/java/javaclasspath.m, scripts/java/javamem.m, scripts/linear-algebra/linsolve.m, scripts/linear-algebra/qzhess.m, scripts/miscellaneous/debug.m, scripts/miscellaneous/desktop.m, scripts/miscellaneous/dir.m, scripts/miscellaneous/dos.m, scripts/miscellaneous/edit.m, scripts/miscellaneous/fact.m, scripts/miscellaneous/getappdata.m, scripts/miscellaneous/inputname.m, scripts/miscellaneous/license.m, scripts/miscellaneous/ls_command.m, scripts/miscellaneous/run.m, scripts/miscellaneous/setfield.m, scripts/miscellaneous/unix.m, scripts/miscellaneous/ver.m, scripts/mk-pkg-add, scripts/mkdoc.pl, scripts/optimization/fminsearch.m, scripts/optimization/optimset.m, scripts/optimization/sqp.m, scripts/pkg/pkg.m, scripts/pkg/private/create_pkgadddel.m, scripts/pkg/private/fix_depends.m, scripts/pkg/private/install.m, scripts/plot/appearance/axis.m, scripts/plot/appearance/box.m, scripts/plot/appearance/clabel.m, scripts/plot/appearance/daspect.m, scripts/plot/appearance/datetick.m, scripts/plot/appearance/grid.m, scripts/plot/appearance/legend.m, scripts/plot/appearance/orient.m, scripts/plot/appearance/shading.m, scripts/plot/appearance/text.m, scripts/plot/appearance/title.m, scripts/plot/appearance/xlabel.m, scripts/plot/appearance/ylabel.m, scripts/plot/appearance/zlabel.m, scripts/plot/draw/area.m, scripts/plot/draw/bar.m, scripts/plot/draw/barh.m, scripts/plot/draw/colorbar.m, scripts/plot/draw/contour.m, scripts/plot/draw/contour3.m, scripts/plot/draw/contourf.m, scripts/plot/draw/ellipsoid.m, scripts/plot/draw/errorbar.m, scripts/plot/draw/ezcontour.m, scripts/plot/draw/ezcontourf.m, scripts/plot/draw/ezmesh.m, scripts/plot/draw/ezpolar.m, scripts/plot/draw/fill.m, scripts/plot/draw/fplot.m, scripts/plot/draw/hist.m, scripts/plot/draw/meshc.m, scripts/plot/draw/meshz.m, scripts/plot/draw/pareto.m, scripts/plot/draw/patch.m, scripts/plot/draw/peaks.m, scripts/plot/draw/pie.m, scripts/plot/draw/pie3.m, scripts/plot/draw/plot.m, scripts/plot/draw/plotyy.m, scripts/plot/draw/private/__bar__.m, scripts/plot/draw/private/__contour__.m, scripts/plot/draw/private/__errplot__.m, scripts/plot/draw/private/__ezplot__.m, scripts/plot/draw/private/__patch__.m, scripts/plot/draw/private/__stem__.m, scripts/plot/draw/rectangle.m, scripts/plot/draw/ribbon.m, scripts/plot/draw/rose.m, scripts/plot/draw/scatter.m, scripts/plot/draw/scatter3.m, scripts/plot/draw/semilogx.m, scripts/plot/draw/shrinkfaces.m, scripts/plot/draw/sombrero.m, scripts/plot/draw/sphere.m, scripts/plot/draw/stairs.m, scripts/plot/draw/stem.m, scripts/plot/draw/stemleaf.m, scripts/plot/draw/surf.m, scripts/plot/draw/surface.m, scripts/plot/draw/surfc.m, scripts/plot/draw/surfl.m, scripts/plot/draw/surfnorm.m, scripts/plot/draw/tetramesh.m, scripts/plot/draw/trimesh.m, scripts/plot/draw/triplot.m, scripts/plot/draw/trisurf.m, scripts/plot/util/__gnuplot_drawnow__.m, scripts/plot/util/__plt_get_axis_arg__.m, scripts/plot/util/axes.m, scripts/plot/util/clf.m, scripts/plot/util/copyobj.m, scripts/plot/util/figure.m, scripts/plot/util/gcbo.m, scripts/plot/util/graphics_toolkit.m, scripts/plot/util/hggroup.m, scripts/plot/util/meshgrid.m, scripts/plot/util/newplot.m, scripts/plot/util/print.m, scripts/plot/util/private/__add_default_menu__.m, scripts/plot/util/private/__fltk_print__.m, scripts/plot/util/private/__gnuplot_print__.m, scripts/plot/util/private/__print_parse_opts__.m, scripts/plot/util/refreshdata.m, scripts/plot/util/subplot.m, scripts/polynomial/conv.m, scripts/polynomial/poly.m, scripts/polynomial/polyeig.m, scripts/polynomial/polyfit.m, scripts/polynomial/polyval.m, scripts/polynomial/private/__splinefit__.m, scripts/polynomial/spline.m, scripts/prefs/prefdir.m, scripts/prefs/preferences.m, scripts/prefs/private/prefsfile.m, scripts/prefs/rmpref.m, scripts/signal/freqz.m, scripts/signal/module.mk, scripts/sparse/eigs.m, scripts/sparse/pcg.m, scripts/sparse/private/__sprand_impl__.m, scripts/sparse/sprand.m, scripts/sparse/sprandn.m, scripts/sparse/spy.m, scripts/sparse/svds.m, scripts/specfun/expint.m, scripts/specfun/factor.m, scripts/special-matrix/gallery.m, scripts/special-matrix/hankel.m, scripts/special-matrix/toeplitz.m, scripts/startup/inputrc, scripts/statistics/base/kurtosis.m, scripts/statistics/base/moment.m, scripts/statistics/base/qqplot.m, scripts/statistics/base/var.m, scripts/statistics/distributions/betarnd.m, scripts/statistics/distributions/binoinv.m, scripts/statistics/distributions/binopdf.m, scripts/statistics/distributions/binornd.m, scripts/statistics/distributions/cauchy_rnd.m, scripts/statistics/distributions/chi2rnd.m, scripts/statistics/distributions/discrete_pdf.m, scripts/statistics/distributions/discrete_rnd.m, scripts/statistics/distributions/empirical_rnd.m, scripts/statistics/distributions/exprnd.m, scripts/statistics/distributions/frnd.m, scripts/statistics/distributions/gamrnd.m, scripts/statistics/distributions/geornd.m, scripts/statistics/distributions/hygernd.m, scripts/statistics/distributions/kolmogorov_smirnov_cdf.m, scripts/statistics/distributions/laplace_cdf.m, scripts/statistics/distributions/laplace_pdf.m, scripts/statistics/distributions/logistic_cdf.m, scripts/statistics/distributions/logistic_pdf.m, scripts/statistics/distributions/lognrnd.m, scripts/statistics/distributions/nbincdf.m, scripts/statistics/distributions/nbininv.m, scripts/statistics/distributions/nbinpdf.m, scripts/statistics/distributions/nbinrnd.m, scripts/statistics/distributions/normrnd.m, scripts/statistics/distributions/poissinv.m, scripts/statistics/distributions/poissrnd.m, scripts/statistics/distributions/tinv.m, scripts/statistics/distributions/trnd.m, scripts/statistics/distributions/unidcdf.m, scripts/statistics/distributions/unidpdf.m, scripts/statistics/distributions/unidrnd.m, scripts/statistics/distributions/unifrnd.m, scripts/statistics/distributions/wblrnd.m, scripts/statistics/models/module.mk, scripts/statistics/tests/kruskal_wallis_test.m, scripts/strings/base2dec.m, scripts/strings/deblank.m, scripts/strings/dec2base.m, scripts/strings/dec2bin.m, scripts/strings/dec2hex.m, scripts/strings/mat2str.m, scripts/strings/ostrsplit.m, scripts/strings/regexptranslate.m, scripts/strings/str2num.m, scripts/strings/strcat.m, scripts/strings/strjoin.m, scripts/strings/strsplit.m, scripts/strings/strtok.m, scripts/strings/strtrim.m, scripts/strings/strtrunc.m, scripts/strings/substr.m, scripts/testfun/__run_test_suite__.m, scripts/testfun/speed.m, scripts/testfun/test.m, scripts/time/asctime.m, scripts/time/datenum.m, scripts/time/datevec.m, scripts/time/weekday.m, src/Makefile.am, test/Makefile.am, test/build-bc-overload-tests.sh, test/build-sparse-tests.sh, test/jit.tst, test/line-continue.tst: Strip trailing whitespace.
author John W. Eaton <jwe@octave.org>
date Tue, 20 Jan 2015 08:26:57 -0500
parents 322eb69e30ad
children 0e1f5a750d00
comparison
equal deleted inserted replaced
19488:8dbd55742112 19627:446c46af4b42
4 @c 4 @c
5 @c Octave is free software; you can redistribute it and/or modify it 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 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 7 @c Free Software Foundation; either version 3 of the License, or (at
8 @c your option) any later version. 8 @c your option) any later version.
9 @c 9 @c
10 @c Octave is distributed in the hope that it will be useful, but WITHOUT 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 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 12 @c FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13 @c for more details. 13 @c for more details.
14 @c 14 @c
15 @c You should have received a copy of the GNU General Public License 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 16 @c along with Octave; see the file COPYING. If not, see
17 @c <http://www.gnu.org/licenses/>. 17 @c <http://www.gnu.org/licenses/>.
18 18
19 @ifhtml 19 @ifhtml
37 @section Creation and Manipulation of Sparse Matrices 37 @section Creation and Manipulation of Sparse Matrices
38 38
39 The size of mathematical problems that can be treated at any particular 39 The size of mathematical problems that can be treated at any particular
40 time is generally limited by the available computing resources. Both, 40 time is generally limited by the available computing resources. Both,
41 the speed of the computer and its available memory place limitation on 41 the speed of the computer and its available memory place limitation on
42 the problem size. 42 the problem size.
43 43
44 There are many classes of mathematical problems which give rise to 44 There are many classes of mathematical problems which give rise to
45 matrices, where a large number of the elements are zero. In this case 45 matrices, where a large number of the elements are zero. In this case
46 it makes sense to have a special matrix type to handle this class of 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 47 problems where only the non-zero elements of the matrix are
66 @subsection Storage of Sparse Matrices 66 @subsection Storage of Sparse Matrices
67 67
68 It is not strictly speaking necessary for the user to understand how 68 It is not strictly speaking necessary for the user to understand how
69 sparse matrices are stored. However, such an understanding will help 69 sparse matrices are stored. However, such an understanding will help
70 to get an understanding of the size of sparse matrices. Understanding 70 to get an understanding of the size of sparse matrices. Understanding
71 the storage technique is also necessary for those users wishing to 71 the storage technique is also necessary for those users wishing to
72 create their own oct-files. 72 create their own oct-files.
73 73
74 There are many different means of storing sparse matrix data. What all 74 There are many different means of storing sparse matrix data. What all
75 of the methods have in common is that they attempt to reduce the complexity 75 of the methods have in common is that they attempt to reduce the complexity
76 and storage given a priori knowledge of the particular class of problems 76 and storage given a priori knowledge of the particular class of problems
77 that will be solved. A good summary of the available techniques for storing 77 that will be solved. A good summary of the available techniques for storing
78 sparse matrix is given by Saad @footnote{Y. Saad "SPARSKIT: A basic toolkit 78 sparse matrix is given by Saad @footnote{Y. Saad "SPARSKIT: A basic toolkit
79 for sparse matrix computation", 1994, 79 for sparse matrix computation", 1994,
80 @url{http://www-users.cs.umn.edu/~saad/software/SPARSKIT/paper.ps}}. 80 @url{http://www-users.cs.umn.edu/~saad/software/SPARSKIT/paper.ps}}.
81 With full matrices, knowledge of the point of an element of the matrix 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. 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 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. 84 of the non-zero elements of the matrix must equally be stored.
85 85
86 An obvious way to do this is by storing the elements of the matrix as 86 An obvious way to do this is by storing the elements of the matrix as
87 triplets, with two elements being their position in the array 87 triplets, with two elements being their position in the array
88 (rows and column) and the third being the data itself. This is conceptually 88 (rows and column) and the third being the data itself. This is conceptually
89 easy to grasp, but requires more storage than is strictly needed. 89 easy to grasp, but requires more storage than is strictly needed.
90 90
91 The storage technique used within Octave is the compressed column 91 The storage technique used within Octave is the compressed column
92 format. It is similar to the Yale format. 92 format. It is similar to the Yale format.
107 107
108 @example 108 @example
109 @group 109 @group
110 for (j = 0; j < nc; j++) 110 for (j = 0; j < nc; j++)
111 for (i = cidx(j); i < cidx(j+1); i++) 111 for (i = cidx(j); i < cidx(j+1); i++)
112 printf ("non-zero element (%i,%i) is %d\n", 112 printf ("non-zero element (%i,%i) is %d\n",
113 ridx(i), j, data(i)); 113 ridx(i), j, data(i));
114 @end group 114 @end group
115 @end example 115 @end example
116 116
117 A clear understanding might be had by considering an example of how the 117 A clear understanding might be had by considering an example of how the
147 @var{data} = [1, 2, 3, 4] 147 @var{data} = [1, 2, 3, 4]
148 @end group 148 @end group
149 @end example 149 @end example
150 150
151 Note that this is the representation of these elements with the first row 151 Note that this is the representation of these elements with the first row
152 and column assumed to start at zero, while in Octave itself the row and 152 and column assumed to start at zero, while in Octave itself the row and
153 column indexing starts at one. Thus the number of elements in the 153 column 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) - 154 @var{i}-th column is given by @code{@var{cidx} (@var{i} + 1) -
155 @var{cidx} (@var{i})}. 155 @var{cidx} (@var{i})}.
156 156
157 Although Octave uses a compressed column format, it should be noted 157 Although Octave uses a compressed column format, it should be noted
158 that compressed row formats are equally possible. However, in the 158 that compressed row formats are equally possible. However, in the
159 context of mixed operations between mixed sparse and dense matrices, 159 context of mixed operations between mixed sparse and dense matrices,
160 it makes sense that the elements of the sparse matrices are in the 160 it makes sense that the elements of the sparse matrices are in the
161 same order as the dense matrices. Octave stores dense matrices in 161 same order as the dense matrices. Octave stores dense matrices in
162 column major ordering, and so sparse matrices are equally stored in 162 column major ordering, and so sparse matrices are equally stored in
163 this manner. 163 this manner.
164 164
165 A further constraint on the sparse matrix storage used by Octave is that 165 A further constraint on the sparse matrix storage used by Octave is that
166 all elements in the rows are stored in increasing order of their row 166 all elements in the rows are stored in increasing order of their row
167 index, which makes certain operations faster. However, it imposes 167 index, which makes certain operations faster. However, it imposes
168 the need to sort the elements on the creation of sparse matrices. Having 168 the need to sort the elements on the creation of sparse matrices. Having
169 disordered elements is potentially an advantage in that it makes operations 169 disordered elements is potentially an advantage in that it makes operations
170 such as concatenating two sparse matrices together easier and faster, however 170 such as concatenating two sparse matrices together easier and faster, however
179 @item Returned from a function 179 @item Returned from a function
180 There are many functions that directly return sparse matrices. These include 180 There are many functions that directly return sparse matrices. These include
181 @dfn{speye}, @dfn{sprand}, @dfn{diag}, etc. 181 @dfn{speye}, @dfn{sprand}, @dfn{diag}, etc.
182 182
183 @item Constructed from matrices or vectors 183 @item Constructed from matrices or vectors
184 The function @dfn{sparse} allows a sparse matrix to be constructed from 184 The function @dfn{sparse} allows a sparse matrix to be constructed from
185 three vectors representing the row, column and data. Alternatively, the 185 three vectors representing the row, column and data. Alternatively, the
186 function @dfn{spconvert} uses a three column matrix format to allow easy 186 function @dfn{spconvert} uses a three column matrix format to allow easy
187 importation of data from elsewhere. 187 importation of data from elsewhere.
188 188
189 @item Created and then filled 189 @item Created and then filled
208 creates an @var{r}-by-@var{c} sparse matrix with a density of filled 208 creates an @var{r}-by-@var{c} sparse matrix with a density of filled
209 elements of @var{d}. 209 elements of @var{d}.
210 210
211 Other functions of interest that directly create sparse matrices, are 211 Other functions of interest that directly create sparse matrices, are
212 @dfn{diag} or its generalization @dfn{spdiags}, that can take the 212 @dfn{diag} or its generalization @dfn{spdiags}, that can take the
213 definition of the diagonals of the matrix and create the sparse matrix 213 definition of the diagonals of the matrix and create the sparse matrix
214 that corresponds to this. For example, 214 that corresponds to this. For example,
215 215
216 @example 216 @example
217 s = diag (sparse (randn (1,n)), -1); 217 s = diag (sparse (randn (1,n)), -1);
218 @end example 218 @end example
231 231
232 @DOCSTRING(sprandn) 232 @DOCSTRING(sprandn)
233 233
234 @DOCSTRING(sprandsym) 234 @DOCSTRING(sprandsym)
235 235
236 The recommended way for the user to create a sparse matrix, is to create 236 The recommended way for the user to create a sparse matrix, is to create
237 two vectors containing the row and column index of the data and a third 237 two vectors containing the row and column index of the data and a third
238 vector of the same size containing the data to be stored. For example, 238 vector of the same size containing the data to be stored. For example,
239 239
240 @example 240 @example
241 @group 241 @group
257 make the creation of the sparse matrix faster. 257 make the creation of the sparse matrix faster.
258 258
259 The function @dfn{spconvert} takes a three or four column real matrix. 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 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 261 the third and four columns, the real and imaginary parts of the sparse
262 matrix. The matrix can contain zero elements and the elements can be 262 matrix. The matrix can contain zero elements and the elements can be
263 sorted in any order. Adding zero elements is a convenient way to define 263 sorted in any order. Adding zero elements is a convenient way to define
264 the size of the sparse matrix. For example: 264 the size of the sparse matrix. For example:
265 265
266 @example 266 @example
267 @group 267 @group
288 @end group 288 @end group
289 @end example 289 @end example
290 290
291 It should be noted, that due to the way that the Octave 291 It should be noted, that due to the way that the Octave
292 assignment functions are written that the assignment will reallocate 292 assignment functions are written that the assignment will reallocate
293 the memory used by the sparse matrix at each iteration of the above loop. 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 294 Therefore the @dfn{spalloc} function ignores the @var{nz} argument and
295 does not pre-assign the memory for the matrix. Therefore, it is vitally 295 does not pre-assign the memory for the matrix. Therefore, it is vitally
296 important that code using to above structure should be vectorized 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 297 as much as possible to minimize the number of assignments and reduce the
298 number of memory allocations. 298 number of memory allocations.
299 299
348 to use of the div or ldiv operators. For example, 348 to use of the div or ldiv operators. For example,
349 349
350 @example 350 @example
351 @group 351 @group
352 a = tril (sprandn (1024, 1024, 0.02), -1) ... 352 a = tril (sprandn (1024, 1024, 0.02), -1) ...
353 + speye (1024); 353 + speye (1024);
354 matrix_type (a); 354 matrix_type (a);
355 ans = Lower 355 ans = Lower
356 @end group 356 @end group
357 @end example 357 @end example
358 358
408 408
409 @noindent 409 @noindent
410 which creates an adjacency matrix @code{A} where node 1 is connected 410 which creates an adjacency matrix @code{A} where node 1 is connected
411 to nodes 2 and 6, node 2 with nodes 1 and 3, etc. The coordinates of 411 to nodes 2 and 6, node 2 with nodes 1 and 3, etc. The coordinates of
412 the nodes are given in the n-by-2 matrix @code{xy}. 412 the nodes are given in the n-by-2 matrix @code{xy}.
413 @ifset htmltex 413 @ifset htmltex
414 @xref{fig:gplot}. 414 @xref{fig:gplot}.
415 415
416 @float Figure,fig:gplot 416 @float Figure,fig:gplot
417 @center @image{gplot,4in} 417 @center @image{gplot,4in}
418 @caption{Simple use of the @dfn{gplot} command.} 418 @caption{Simple use of the @dfn{gplot} command.}
452 452
453 Many Octave functions have been overloaded to work with either sparse or full 453 Many Octave functions have been overloaded to work with either sparse or full
454 matrices. There is no difference in calling convention when using an 454 matrices. There is no difference in calling convention when using an
455 overloaded function with a sparse matrix, however, there is also no access to 455 overloaded function with a sparse matrix, however, there is also no access to
456 potentially sparse-specific features. At any time the sparse matrix specific 456 potentially sparse-specific features. At any time the sparse matrix specific
457 version of a function can be used by explicitly calling its function name. 457 version of a function can be used by explicitly calling its function name.
458 458
459 The table below lists all of the sparse functions of Octave. Note that the 459 The table below lists all of the sparse functions of Octave. Note that the
460 names of the specific sparse forms of the functions are typically the same as 460 names of the specific sparse forms of the functions are typically the same as
461 the general versions with a @dfn{sp} prefix. In the table below, and in the 461 the general versions with a @dfn{sp} prefix. In the table below, and in the
462 rest of this article, the specific sparse versions of functions are used. 462 rest of this article, the specific sparse versions of functions are used.
463 463
464 @c Table includes in comments the missing sparse functions 464 @c Table includes in comments the missing sparse functions
465 465
466 @table @asis 466 @table @asis
467 @item Generate sparse matrices: 467 @item Generate sparse matrices:
468 @dfn{spalloc}, @dfn{spdiags}, @dfn{speye}, @dfn{sprand}, 468 @dfn{spalloc}, @dfn{spdiags}, @dfn{speye}, @dfn{sprand},
469 @dfn{sprandn}, @dfn{sprandsym} 469 @dfn{sprandn}, @dfn{sprandsym}
470 470
471 @item Sparse matrix conversion: 471 @item Sparse matrix conversion:
472 @dfn{full}, @dfn{sparse}, @dfn{spconvert} 472 @dfn{full}, @dfn{sparse}, @dfn{spconvert}
473 473
474 @item Manipulate sparse matrices 474 @item Manipulate sparse matrices
475 @dfn{issparse}, @dfn{nnz}, @dfn{nonzeros}, @dfn{nzmax}, 475 @dfn{issparse}, @dfn{nnz}, @dfn{nonzeros}, @dfn{nzmax},
476 @dfn{spfun}, @dfn{spones}, @dfn{spy} 476 @dfn{spfun}, @dfn{spones}, @dfn{spy}
477 477
478 @item Graph Theory: 478 @item Graph Theory:
479 @dfn{etree}, @dfn{etreeplot}, @dfn{gplot}, 479 @dfn{etree}, @dfn{etreeplot}, @dfn{gplot},
480 @dfn{treeplot} 480 @dfn{treeplot}
481 @c @dfn{treelayout} 481 @c @dfn{treelayout}
482 482
483 @item Sparse matrix reordering: 483 @item Sparse matrix reordering:
484 @dfn{amd}, @dfn{ccolamd}, @dfn{colamd}, @dfn{colperm}, @dfn{csymamd}, 484 @dfn{amd}, @dfn{ccolamd}, @dfn{colamd}, @dfn{colperm}, @dfn{csymamd},
488 @dfn{condest}, @dfn{eigs}, @dfn{matrix_type}, @dfn{normest}, @dfn{sprank}, 488 @dfn{condest}, @dfn{eigs}, @dfn{matrix_type}, @dfn{normest}, @dfn{sprank},
489 @dfn{spaugment}, @dfn{svds} 489 @dfn{spaugment}, @dfn{svds}
490 490
491 @item Iterative techniques: 491 @item Iterative techniques:
492 @dfn{luinc}, @dfn{pcg}, @dfn{pcr} 492 @dfn{luinc}, @dfn{pcg}, @dfn{pcr}
493 @c @dfn{bicg}, @dfn{bicgstab}, @dfn{cholinc}, @dfn{cgs}, @dfn{gmres}, 493 @c @dfn{bicg}, @dfn{bicgstab}, @dfn{cholinc}, @dfn{cgs}, @dfn{gmres},
494 @c @dfn{lsqr}, @dfn{minres}, @dfn{qmr}, @dfn{symmlq} 494 @c @dfn{lsqr}, @dfn{minres}, @dfn{qmr}, @dfn{symmlq}
495 495
496 @item Miscellaneous: 496 @item Miscellaneous:
497 @dfn{spparms}, @dfn{symbfact}, @dfn{spstats} 497 @dfn{spparms}, @dfn{symbfact}, @dfn{spstats}
498 @end table 498 @end table
504 details. 504 details.
505 505
506 @node Return Types of Operators and Functions 506 @node Return Types of Operators and Functions
507 @subsubsection Return Types of Operators and Functions 507 @subsubsection Return Types of Operators and Functions
508 508
509 The two basic reasons to use sparse matrices are to reduce the memory 509 The two basic reasons to use sparse matrices are to reduce the memory
510 usage and to not have to do calculations on zero elements. The two are 510 usage and to not have to do calculations on zero elements. The two are
511 closely related in that the computation time on a sparse matrix operator 511 closely related in that the computation time on a sparse matrix operator
512 or function is roughly linear with the number of non-zero elements. 512 or function is roughly linear with the number of non-zero elements.
513 513
514 Therefore, there is a certain density of non-zero elements of a matrix 514 Therefore, there is a certain density of non-zero elements of a matrix
515 where it no longer makes sense to store it as a sparse matrix, but rather 515 where it no longer makes sense to store it as a sparse matrix, but rather
516 as a full matrix. For this reason operators and functions that have a 516 as a full matrix. For this reason operators and functions that have a
517 high probability of returning a full matrix will always return one. For 517 high probability of returning a full matrix will always return one. For
518 example adding a scalar constant to a sparse matrix will almost always 518 example adding a scalar constant to a sparse matrix will almost always
519 make it a full matrix, and so the example, 519 make it a full matrix, and so the example,
520 520
521 @example 521 @example
526 0 0 1 526 0 0 1
527 @end group 527 @end group
528 @end example 528 @end example
529 529
530 @noindent 530 @noindent
531 returns a full matrix as can be seen. 531 returns a full matrix as can be seen.
532 532
533 533
534 Additionally, if @code{sparse_auto_mutate} is true, all sparse functions 534 Additionally, if @code{sparse_auto_mutate} is true, all sparse functions
535 test the amount of memory occupied by the sparse matrix to see if the 535 test the amount of memory occupied by the sparse matrix to see if the
536 amount of storage used is larger than the amount used by the full 536 amount of storage used is larger than the amount used by the full
537 equivalent. Therefore @code{speye (2) * 1} will return a full matrix as 537 equivalent. Therefore @code{speye (2) * 1} will return a full matrix as
538 the memory used is smaller for the full version than the sparse version. 538 the memory used is smaller for the full version than the sparse version.
539 539
540 As all of the mixed operators and functions between full and sparse 540 As all of the mixed operators and functions between full and sparse
541 matrices exist, in general this does not cause any problems. However, 541 matrices exist, in general this does not cause any problems. However,
542 one area where it does cause a problem is where a sparse matrix is 542 one area where it does cause a problem is where a sparse matrix is
543 promoted to a full matrix, where subsequent operations would resparsify 543 promoted to a full matrix, where subsequent operations would resparsify
544 the matrix. Such cases are rare, but can be artificially created, for 544 the matrix. Such cases are rare, but can be artificially created, for
545 example @code{(fliplr (speye (3)) + speye (3)) - speye (3)} gives a full 545 example @code{(fliplr (speye (3)) + speye (3)) - speye (3)} gives a full
546 matrix when it should give a sparse one. In general, where such cases 546 matrix when it should give a sparse one. In general, where such cases
547 occur, they impose only a small memory penalty. 547 occur, they impose only a small memory penalty.
548 548
549 There is however one known case where this behavior of Octave's 549 There is however one known case where this behavior of Octave's
550 sparse matrices will cause a problem. That is in the handling of the 550 sparse 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 551 @dfn{diag} function. Whether @dfn{diag} returns a sparse or full matrix
552 depending on the type of its input arguments. So 552 depending on the type of its input arguments. So
553 553
554 @example 554 @example
555 a = diag (sparse ([1,2,3]), -1); 555 a = diag (sparse ([1,2,3]), -1);
556 @end example 556 @end example
557 557
558 @noindent 558 @noindent
559 should return a sparse matrix. To ensure this actually happens, the 559 should return a sparse matrix. To ensure this actually happens, the
560 @dfn{sparse} function, and other functions based on it like @dfn{speye}, 560 @dfn{sparse} function, and other functions based on it like @dfn{speye},
561 always returns a sparse matrix, even if the memory used will be larger 561 always returns a sparse matrix, even if the memory used will be larger
562 than its full representation. 562 than its full representation.
563 563
564 @DOCSTRING(sparse_auto_mutate) 564 @DOCSTRING(sparse_auto_mutate)
565 565
566 Note that the @code{sparse_auto_mutate} option is incompatible with 566 Note that the @code{sparse_auto_mutate} option is incompatible with
571 571
572 The attempt has been made to make sparse matrices behave in exactly the 572 The attempt has been made to make sparse matrices behave in exactly the
573 same manner as there full counterparts. However, there are certain differences 573 same manner as there full counterparts. However, there are certain differences
574 and especially differences with other products sparse implementations. 574 and especially differences with other products sparse implementations.
575 575
576 First, the @qcode{"./"} and @qcode{".^"} operators must be used with care. 576 First, the @qcode{"./"} and @qcode{".^"} operators must be used with care.
577 Consider what the examples 577 Consider what the examples
578 578
579 @example 579 @example
580 @group 580 @group
581 s = speye (4); 581 s = speye (4);
590 590
591 @noindent 591 @noindent
592 will give. The first example of @var{s} raised to the power of 2 causes 592 will give. The first example of @var{s} raised to the power of 2 causes
593 no problems. However @var{s} raised element-wise to itself involves a 593 no problems. However @var{s} raised element-wise to itself involves a
594 large number of terms @code{0 .^ 0} which is 1. There @code{@var{s} .^ 594 large number of terms @code{0 .^ 0} which is 1. There @code{@var{s} .^
595 @var{s}} is a full matrix. 595 @var{s}} is a full matrix.
596 596
597 Likewise @code{@var{s} .^ -2} involves terms like @code{0 .^ -2} which 597 Likewise @code{@var{s} .^ -2} involves terms like @code{0 .^ -2} which
598 is infinity, and so @code{@var{s} .^ -2} is equally a full matrix. 598 is infinity, and so @code{@var{s} .^ -2} is equally a full matrix.
599 599
600 For the "./" operator @code{@var{s} ./ 2} has no problems, but 600 For the "./" operator @code{@var{s} ./ 2} has no problems, but
601 @code{2 ./ @var{s}} involves a large number of infinity terms as well 601 @code{2 ./ @var{s}} involves a large number of infinity terms as well
602 and is equally a full matrix. The case of @code{@var{s} ./ @var{s}} 602 and is equally a full matrix. The case of @code{@var{s} ./ @var{s}}
603 involves terms like @code{0 ./ 0} which is a @code{NaN} and so this 603 involves terms like @code{0 ./ 0} which is a @code{NaN} and so this
604 is equally a full matrix with the zero elements of @var{s} filled with 604 is equally a full matrix with the zero elements of @var{s} filled with
605 @code{NaN} values. 605 @code{NaN} values.
606 606
607 The above behavior is consistent with full matrices, but is not 607 The above behavior is consistent with full matrices, but is not
608 consistent with sparse implementations in other products. 608 consistent with sparse implementations in other products.
609 609
610 A particular problem of sparse matrices comes about due to the fact that 610 A particular problem of sparse matrices comes about due to the fact that
611 as the zeros are not stored, the sign-bit of these zeros is equally not 611 as the zeros are not stored, the sign-bit of these zeros is equally not
612 stored. In certain cases the sign-bit of zero is important. For example: 612 stored. In certain cases the sign-bit of zero is important. For example:
620 c = 1 ./ sparse (a) 620 c = 1 ./ sparse (a)
621 @result{} Inf Inf 621 @result{} Inf Inf
622 Inf Inf 622 Inf Inf
623 @end group 623 @end group
624 @end example 624 @end example
625 625
626 To correct this behavior would mean that zero elements with a negative 626 To correct this behavior would mean that zero elements with a negative
627 sign-bit would need to be stored in the matrix to ensure that their 627 sign-bit would need to be stored in the matrix to ensure that their
628 sign-bit was respected. This is not done at this time, for reasons of 628 sign-bit was respected. This is not done at this time, for reasons of
629 efficiency, and so the user is warned that calculations where the sign-bit 629 efficiency, and so the user is warned that calculations where the sign-bit
630 of zero is important must not be done using sparse matrices. 630 of zero is important must not be done using sparse matrices.
631 631
632 In general any function or operator used on a sparse matrix will 632 In general any function or operator used on a sparse matrix will
643 then @dfn{symamd} or @dfn{csymamd} should be used. Otherwise 643 then @dfn{symamd} or @dfn{csymamd} should be used. Otherwise
644 @dfn{amd}, @dfn{colamd} or @dfn{ccolamd} should be used. For completeness 644 @dfn{amd}, @dfn{colamd} or @dfn{ccolamd} should be used. For completeness
645 the reordering functions @dfn{colperm} and @dfn{randperm} are 645 the reordering functions @dfn{colperm} and @dfn{randperm} are
646 also available. 646 also available.
647 647
648 @xref{fig:simplematrix}, for an example of the structure of a simple 648 @xref{fig:simplematrix}, for an example of the structure of a simple
649 positive definite matrix. 649 positive definite matrix.
650 650
651 @float Figure,fig:simplematrix 651 @float Figure,fig:simplematrix
652 @center @image{spmatrix,4in} 652 @center @image{spmatrix,4in}
653 @caption{Structure of simple sparse matrix.} 653 @caption{Structure of simple sparse matrix.}
654 @end float 654 @end float
655 655
656 The standard Cholesky@tie{}factorization of this matrix can be 656 The standard Cholesky@tie{}factorization of this matrix can be
657 obtained by the same command that would be used for a full 657 obtained by the same command that would be used for a full
658 matrix. This can be visualized with the command 658 matrix. This can be visualized with the command
659 @code{r = chol (A); spy (r);}. 659 @code{r = chol (A); spy (r);}.
660 @xref{fig:simplechol}. 660 @xref{fig:simplechol}.
661 The original matrix had 661 The original matrix had
662 @ifinfo 662 @ifinfo
663 @ifnothtml 663 @ifnothtml
664 43 664 43
665 @end ifnothtml 665 @end ifnothtml
666 @end ifinfo 666 @end ifinfo
676 @ifset htmltex 676 @ifset htmltex
677 10200, 677 10200,
678 @end ifset 678 @end ifset
679 with only half of the symmetric matrix being stored. This 679 with only half of the symmetric matrix being stored. This
680 is a significant level of fill in, and although not an issue 680 is a significant level of fill in, and although not an issue
681 for such a small test case, can represents a large overhead 681 for such a small test case, can represents a large overhead
682 in working with other sparse matrices. 682 in working with other sparse matrices.
683 683
684 The appropriate sparsity preserving permutation of the original 684 The appropriate sparsity preserving permutation of the original
685 matrix is given by @dfn{symamd} and the factorization using this 685 matrix is given by @dfn{symamd} and the factorization using this
686 reordering can be visualized using the command @code{q = symamd (A); 686 reordering can be visualized using the command @code{q = symamd (A);
687 r = chol (A(q,q)); spy (r)}. This gives 687 r = chol (A(q,q)); spy (r)}. This gives
688 @ifinfo 688 @ifinfo
689 @ifnothtml 689 @ifnothtml
690 29 690 29
691 @end ifnothtml 691 @end ifnothtml
692 @end ifinfo 692 @end ifinfo
736 @DOCSTRING(symrcm) 736 @DOCSTRING(symrcm)
737 737
738 @node Sparse Linear Algebra 738 @node Sparse Linear Algebra
739 @section Linear Algebra on Sparse Matrices 739 @section Linear Algebra on Sparse Matrices
740 740
741 Octave includes a polymorphic solver for sparse matrices, where 741 Octave includes a polymorphic solver for sparse matrices, where
742 the exact solver used to factorize the matrix, depends on the properties 742 the exact solver used to factorize the matrix, depends on the properties
743 of the sparse matrix itself. Generally, the cost of determining the matrix type 743 of the sparse matrix itself. Generally, the cost of determining the matrix type
744 is small relative to the cost of factorizing the matrix itself, but in any 744 is small relative to the cost of factorizing the matrix itself, but in any
745 case the matrix type is cached once it is calculated, so that it is not 745 case the matrix type is cached once it is calculated, so that it is not
746 re-determined each time it is used in a linear equation. 746 re-determined each time it is used in a linear equation.
755 755
756 @item If the matrix is square, banded and if the band density is less 756 @item If the matrix is square, banded and if the band density is less
757 than that given by @code{spparms ("bandden")} continue, else goto 4. 757 than that given by @code{spparms ("bandden")} continue, else goto 4.
758 758
759 @enumerate a 759 @enumerate a
760 @item If the matrix is tridiagonal and the right-hand side is not sparse 760 @item If the matrix is tridiagonal and the right-hand side is not sparse
761 continue, else goto 3b. 761 continue, else goto 3b.
762 762
763 @enumerate 763 @enumerate
764 @item If the matrix is Hermitian, with a positive real diagonal, attempt 764 @item If the matrix is Hermitian, with a positive real diagonal, attempt
765 Cholesky@tie{}factorization using @sc{lapack} xPTSV. 765 Cholesky@tie{}factorization using @sc{lapack} xPTSV.
766 766
767 @item If the above failed or the matrix is not Hermitian with a positive 767 @item If the above failed or the matrix is not Hermitian with a positive
768 real diagonal use Gaussian elimination with pivoting using 768 real diagonal use Gaussian elimination with pivoting using
769 @sc{lapack} xGTSV, and goto 8. 769 @sc{lapack} xGTSV, and goto 8.
770 @end enumerate 770 @end enumerate
771 771
772 @item If the matrix is Hermitian with a positive real diagonal, attempt 772 @item If the matrix is Hermitian with a positive real diagonal, attempt
773 Cholesky@tie{}factorization using @sc{lapack} xPBTRF. 773 Cholesky@tie{}factorization using @sc{lapack} xPBTRF.
774 774
775 @item if the above failed or the matrix is not Hermitian with a positive 775 @item if the above failed or the matrix is not Hermitian with a positive
776 real diagonal use Gaussian elimination with pivoting using 776 real diagonal use Gaussian elimination with pivoting using
777 @sc{lapack} xGBTRF, and goto 8. 777 @sc{lapack} xGBTRF, and goto 8.
778 @end enumerate 778 @end enumerate
779 779
780 @item If the matrix is upper or lower triangular perform a sparse forward 780 @item If the matrix is upper or lower triangular perform a sparse forward
781 or backward substitution, and goto 8 781 or backward substitution, and goto 8
782 782
783 @item If the matrix is an upper triangular matrix with column permutations 783 @item If the matrix is an upper triangular matrix with column permutations
784 or lower triangular matrix with row permutations, perform a sparse forward 784 or lower triangular matrix with row permutations, perform a sparse forward
785 or backward substitution, and goto 8 785 or backward substitution, and goto 8
786 786
787 @item If the matrix is square, Hermitian with a real positive diagonal, attempt 787 @item If the matrix is square, Hermitian with a real positive diagonal, attempt
788 sparse Cholesky@tie{}factorization using @sc{cholmod}. 788 sparse Cholesky@tie{}factorization using @sc{cholmod}.
789 789
790 @item If the sparse Cholesky@tie{}factorization failed or the matrix is not 790 @item If the sparse Cholesky@tie{}factorization failed or the matrix is not
791 Hermitian with a real positive diagonal, and the matrix is square, factorize 791 Hermitian with a real positive diagonal, and the matrix is square, factorize
792 using @sc{umfpack}. 792 using @sc{umfpack}.
793 793
794 @item If the matrix is not square, or any of the previous solvers flags 794 @item If the matrix is not square, or any of the previous solvers flags
795 a singular or near singular matrix, find a minimum norm solution using 795 a singular or near singular matrix, find a minimum norm solution using
796 @sc{cxsparse}@footnote{The @sc{cholmod}, @sc{umfpack} and @sc{cxsparse} packages were 796 @sc{cxsparse}@footnote{The @sc{cholmod}, @sc{umfpack} and @sc{cxsparse} packages were
883 partial differential equations that do not have closed form solutions, 883 partial differential equations that do not have closed form solutions,
884 typically because of the complex shape of the domain. 884 typically because of the complex shape of the domain.
885 885
886 In order to motivate this application, we consider the boundary value 886 In order to motivate this application, we consider the boundary value
887 Laplace equation. This system can model scalar potential fields, such 887 Laplace equation. This system can model scalar potential fields, such
888 as heat or electrical potential. Given a medium 888 as heat or electrical potential. Given a medium
889 @tex 889 @tex
890 $\Omega$ with boundary $\partial\Omega$. At all points on the $\partial\Omega$ 890 $\Omega$ with boundary $\partial\Omega$. At all points on the $\partial\Omega$
891 the boundary conditions are known, and we wish to calculate the potential in 891 the boundary conditions are known, and we wish to calculate the potential in
892 $\Omega$. 892 $\Omega$.
893 @end tex 893 @end tex
901 (Neumann boundary condition), or a weighted sum of the potential and 901 (Neumann boundary condition), or a weighted sum of the potential and
902 its derivative (Cauchy boundary condition). 902 its derivative (Cauchy boundary condition).
903 903
904 In a thermal model, we want to calculate the temperature in 904 In a thermal model, we want to calculate the temperature in
905 @tex 905 @tex
906 $\Omega$ 906 $\Omega$
907 @end tex 907 @end tex
908 @ifnottex 908 @ifnottex
909 Omega 909 Omega
910 @end ifnottex 910 @end ifnottex
911 and know the boundary temperature (Dirichlet condition) 911 and know the boundary temperature (Dirichlet condition)
912 or heat flux (from which we can calculate the Neumann condition 912 or heat flux (from which we can calculate the Neumann condition
913 by dividing by the thermal conductivity at the boundary). Similarly, 913 by dividing by the thermal conductivity at the boundary). Similarly,
914 in an electrical model, we want to calculate the voltage in 914 in an electrical model, we want to calculate the voltage in
915 @tex 915 @tex
916 $\Omega$ 916 $\Omega$
917 @end tex 917 @end tex
918 @ifnottex 918 @ifnottex
919 Omega 919 Omega
920 @end ifnottex 920 @end ifnottex
921 and know the boundary voltage (Dirichlet) or current 921 and know the boundary voltage (Dirichlet) or current
922 (Neumann condition after diving by the electrical conductivity). 922 (Neumann condition after diving by the electrical conductivity).
923 In an electrical model, it is common for much of the boundary 923 In an electrical model, it is common for much of the boundary
924 to be electrically isolated; this is a Neumann boundary condition 924 to be electrically isolated; this is a Neumann boundary condition
925 with the current equal to zero. 925 with the current equal to zero.
926 926
927 The simplest finite element models will divide 927 The simplest finite element models will divide
928 @tex 928 @tex
929 $\Omega$ 929 $\Omega$
930 @end tex 930 @end tex
931 @ifnottex 931 @ifnottex
932 Omega 932 Omega
933 @end ifnottex 933 @end ifnottex
934 into simplexes (triangles in 2D, pyramids in 3D). 934 into simplexes (triangles in 2D, pyramids in 3D).
935 @ifset htmltex 935 @ifset htmltex
936 We take as a 3-D example a cylindrical liquid filled tank with a small 936 We take as a 3-D example a cylindrical liquid filled tank with a small
937 non-conductive ball from the EIDORS project@footnote{EIDORS - Electrical 937 non-conductive ball from the EIDORS project@footnote{EIDORS - Electrical
938 Impedance Tomography and Diffuse optical Tomography Reconstruction Software 938 Impedance Tomography and Diffuse optical Tomography Reconstruction Software
939 @url{http://eidors3d.sourceforge.net}}. This is model is designed to reflect 939 @url{http://eidors3d.sourceforge.net}}. This is model is designed to reflect
940 an application of electrical impedance tomography, where current patterns 940 an application of electrical impedance tomography, where current patterns
941 are applied to such a tank in order to image the internal conductivity 941 are applied to such a tank in order to image the internal conductivity
942 distribution. In order to describe the FEM geometry, we have a matrix of 942 distribution. In order to describe the FEM geometry, we have a matrix of
943 vertices @code{nodes} and simplices @code{elems}. 943 vertices @code{nodes} and simplices @code{elems}.
944 @end ifset 944 @end ifset
945 945
946 The following example creates a simple rectangular 2-D electrically 946 The following example creates a simple rectangular 2-D electrically
947 conductive medium with 10 V and 20 V imposed on opposite sides 947 conductive medium with 10 V and 20 V imposed on opposite sides
948 (Dirichlet boundary conditions). All other edges are electrically 948 (Dirichlet boundary conditions). All other edges are electrically
949 isolated. 949 isolated.
950 950
951 @example 951 @example
952 @group 952 @group
959 elems = []; 959 elems = [];
960 for idx = 1:w-1 960 for idx = 1:w-1
961 widx = (idx-1)*h; 961 widx = (idx-1)*h;
962 elems = [elems; ... 962 elems = [elems; ...
963 widx+[(1:h-1);(2:h);h+(1:h-1)]'; ... 963 widx+[(1:h-1);(2:h);h+(1:h-1)]'; ...
964 widx+[(2:h);h+(2:h);h+(1:h-1)]' ]; 964 widx+[(2:h);h+(2:h);h+(1:h-1)]' ];
965 endfor 965 endfor
966 966
967 E = size (elems,1); # No. of simplices 967 E = size (elems,1); # No. of simplices
968 N = size (nodes,1); # No. of vertices 968 N = size (nodes,1); # No. of vertices
969 D = size (elems,2); # dimensions+1 969 D = size (elems,2); # dimensions+1
984 2 3 4 5 7 8 9 @dots{} 984 2 3 4 5 7 8 9 @dots{}
985 6 7 8 9 6 7 8 @dots{} 985 6 7 8 9 6 7 8 @dots{}
986 @end group 986 @end group
987 @end example 987 @end example
988 988
989 Using a first order FEM, we approximate the electrical conductivity 989 Using a first order FEM, we approximate the electrical conductivity
990 distribution in 990 distribution in
991 @tex 991 @tex
992 $\Omega$ 992 $\Omega$
993 @end tex 993 @end tex
994 @ifnottex 994 @ifnottex
995 Omega 995 Omega
996 @end ifnottex 996 @end ifnottex
997 as constant on each simplex (represented by the vector @code{conductivity}). 997 as constant on each simplex (represented by the vector @code{conductivity}).
998 Based on the finite element geometry, we first calculate a system (or 998 Based on the finite element geometry, we first calculate a system (or
999 stiffness) matrix for each simplex (represented as 3-by-3 elements on the 999 stiffness) matrix for each simplex (represented as 3-by-3 elements on the
1000 diagonal of the element-wise system matrix @code{SE}. Based on @code{SE} 1000 diagonal of the element-wise system matrix @code{SE}. Based on @code{SE}
1001 and a N-by-DE connectivity matrix @code{C}, representing the connections 1001 and a N-by-DE connectivity matrix @code{C}, representing the connections
1002 between simplices and vertices, the global connectivity matrix @code{S} is 1002 between simplices and vertices, the global connectivity matrix @code{S} is
1003 calculated. 1003 calculated.
1004 1004
1005 @example 1005 @example
1006 ## Element conductivity 1006 ## Element conductivity
1016 ones(1,D) + ones(D*E,1)*(1:D) ; 1016 ones(1,D) + ones(D*E,1)*(1:D) ;
1017 Sjidx = [1:D*E]'*ones (1,D); 1017 Sjidx = [1:D*E]'*ones (1,D);
1018 Sdata = zeros (D*E,D); 1018 Sdata = zeros (D*E,D);
1019 dfact = factorial (D-1); 1019 dfact = factorial (D-1);
1020 for j = 1:E 1020 for j = 1:E
1021 a = inv ([ones(D,1), ... 1021 a = inv ([ones(D,1), ...
1022 nodes(elems(j,:), :)]); 1022 nodes(elems(j,:), :)]);
1023 const = conductivity(j) * 2 / ... 1023 const = conductivity(j) * 2 / ...
1024 dfact / abs (det (a)); 1024 dfact / abs (det (a));
1025 Sdata(D*(j-1)+(1:D),:) = const * ... 1025 Sdata(D*(j-1)+(1:D),:) = const * ...
1026 a(2:D,:)' * a(2:D,:); 1026 a(2:D,:)' * a(2:D,:);
1029 SE = sparse(Siidx,Sjidx,Sdata); 1029 SE = sparse(Siidx,Sjidx,Sdata);
1030 ## Global system matrix 1030 ## Global system matrix
1031 S = C'* SE *C; 1031 S = C'* SE *C;
1032 @end example 1032 @end example
1033 1033
1034 The system matrix acts like the conductivity 1034 The system matrix acts like the conductivity
1035 @tex 1035 @tex
1036 $S$ 1036 $S$
1037 @end tex 1037 @end tex
1038 @ifnottex 1038 @ifnottex
1039 @code{S} 1039 @code{S}
1040 @end ifnottex 1040 @end ifnottex
1041 in Ohm's law 1041 in Ohm's law
1042 @tex 1042 @tex
1043 $SV = I$. 1043 $SV = I$.
1044 @end tex 1044 @end tex
1045 @ifnottex 1045 @ifnottex
1046 @code{S * V = I}. 1046 @code{S * V = I}.
1047 @end ifnottex 1047 @end ifnottex
1048 Based on the Dirichlet and Neumann boundary conditions, we are able to 1048 Based on the Dirichlet and Neumann boundary conditions, we are able to
1049 solve for the voltages at each vertex @code{V}. 1049 solve for the voltages at each vertex @code{V}.
1050 1050
1051 @example 1051 @example
1052 ## Dirichlet boundary conditions 1052 ## Dirichlet boundary conditions
1053 D_nodes = [1:5, 51:55]; 1053 D_nodes = [1:5, 51:55];
1054 D_value = [10*ones(1,5), 20*ones(1,5)]; 1054 D_value = [10*ones(1,5), 20*ones(1,5)];
1055 1055
1056 V = zeros (N,1); 1056 V = zeros (N,1);
1057 V(D_nodes) = D_value; 1057 V(D_nodes) = D_value;
1058 idx = 1:N; # vertices without Dirichlet 1058 idx = 1:N; # vertices without Dirichlet
1059 # boundary condns 1059 # boundary condns
1060 idx(D_nodes) = []; 1060 idx(D_nodes) = [];
1061 1061
1062 ## Neumann boundary conditions. Note that 1062 ## Neumann boundary conditions. Note that
1063 ## N_value must be normalized by the 1063 ## N_value must be normalized by the
1070 1070
1071 V(idx) = S(idx,idx) \ ( Q(idx) - ... 1071 V(idx) = S(idx,idx) \ ( Q(idx) - ...
1072 S(idx,D_nodes) * V(D_nodes)); 1072 S(idx,D_nodes) * V(D_nodes));
1073 @end example 1073 @end example
1074 1074
1075 Finally, in order to display the solution, we show each solved voltage 1075 Finally, in order to display the solution, we show each solved voltage
1076 value in the z-axis for each simplex vertex. 1076 value in the z-axis for each simplex vertex.
1077 @ifset htmltex 1077 @ifset htmltex
1078 @xref{fig:femmodel}. 1078 @xref{fig:femmodel}.
1079 @end ifset 1079 @end ifset
1080 1080
1082 @group 1082 @group
1083 elemx = elems(:,[1,2,3,1])'; 1083 elemx = elems(:,[1,2,3,1])';
1084 xelems = reshape (nodes(elemx, 1), 4, E); 1084 xelems = reshape (nodes(elemx, 1), 4, E);
1085 yelems = reshape (nodes(elemx, 2), 4, E); 1085 yelems = reshape (nodes(elemx, 2), 4, E);
1086 velems = reshape (V(elemx), 4, E); 1086 velems = reshape (V(elemx), 4, E);
1087 plot3 (xelems,yelems,velems,"k"); 1087 plot3 (xelems,yelems,velems,"k");
1088 print "grid.eps"; 1088 print "grid.eps";
1089 @end group 1089 @end group
1090 @end example 1090 @end example
1091 1091
1092 1092
1093 @ifset htmltex 1093 @ifset htmltex
1094 @float Figure,fig:femmodel 1094 @float Figure,fig:femmodel
1095 @center @image{grid,4in} 1095 @center @image{grid,4in}
1096 @caption{Example finite element model the showing triangular elements. 1096 @caption{Example finite element model the showing triangular elements.
1097 The height of each vertex corresponds to the solution value.} 1097 The height of each vertex corresponds to the solution value.}
1098 @end float 1098 @end float
1099 @end ifset 1099 @end ifset