Mercurial > octave
view liboctave/array/CSparse.cc @ 22755:3a2b891d0b33
maint: Standardize Copyright formatting.
* Makefile.am, README, build-aux/check-subst-vars.in.sh,
build-aux/find-files-with-tests.sh, build-aux/mk-builtins.sh,
build-aux/mk-default-qt-settings.in.sh, build-aux/mk-f77-def.in.sh,
build-aux/mk-hg-id.sh, build-aux/mk-mxarray-h.in.sh,
build-aux/mk-octave-config-h.sh, build-aux/mk-opts.pl,
build-aux/mk-version-h.in.sh, build-aux/subst-config-vals.in.sh,
build-aux/subst-cross-config-vals.in.sh, build-aux/subst-default-vals.in.sh,
build-aux/subst-f77-isnan-macro.in.sh, build-aux/subst-script-vals.in.sh,
configure.ac, 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/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/genpropdoc.m,
doc/interpreter/geometry.txi, doc/interpreter/geometryimages.m,
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/linalg.txi,
doc/interpreter/macros.texi, doc/interpreter/matrix.txi,
doc/interpreter/mk-doc-cache.pl, doc/interpreter/mkoctfile.1,
doc/interpreter/nonlin.txi, doc/interpreter/numbers.txi,
doc/interpreter/obsolete.txi, doc/interpreter/octave-cli.1,
doc/interpreter/octave-config.1, doc/interpreter/octave.1,
doc/interpreter/octave.css, 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/plotimages.m, 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/var.txi,
doc/interpreter/vectorize.txi, 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/refcard-a4.tex,
doc/refcard/refcard-legal.tex, doc/refcard/refcard-letter.tex,
doc/refcard/refcard.tex, etc/HACKING, etc/icons/octave.appdata.xml.in,
libgui/graphics/Backend.cc, libgui/graphics/Backend.h,
libgui/graphics/BaseControl.cc, libgui/graphics/BaseControl.h,
libgui/graphics/ButtonControl.cc, libgui/graphics/ButtonControl.h,
libgui/graphics/ButtonGroup.cc, libgui/graphics/ButtonGroup.h,
libgui/graphics/Canvas.cc, libgui/graphics/Canvas.h,
libgui/graphics/CheckBoxControl.cc, libgui/graphics/CheckBoxControl.h,
libgui/graphics/Container.cc, libgui/graphics/Container.h,
libgui/graphics/ContextMenu.cc, libgui/graphics/ContextMenu.h,
libgui/graphics/EditControl.cc, libgui/graphics/EditControl.h,
libgui/graphics/Figure.cc, libgui/graphics/Figure.h,
libgui/graphics/FigureWindow.cc, libgui/graphics/FigureWindow.h,
libgui/graphics/GLCanvas.cc, libgui/graphics/GLCanvas.h,
libgui/graphics/GenericEventNotify.h, libgui/graphics/KeyMap.cc,
libgui/graphics/KeyMap.h, libgui/graphics/ListBoxControl.cc,
libgui/graphics/ListBoxControl.h, libgui/graphics/Logger.cc,
libgui/graphics/Logger.h, libgui/graphics/Menu.cc, libgui/graphics/Menu.h,
libgui/graphics/MenuContainer.h, libgui/graphics/MouseModeActionGroup.cc,
libgui/graphics/MouseModeActionGroup.h, libgui/graphics/Object.cc,
libgui/graphics/Object.h, libgui/graphics/ObjectFactory.cc,
libgui/graphics/ObjectFactory.h, libgui/graphics/ObjectProxy.cc,
libgui/graphics/ObjectProxy.h, libgui/graphics/Panel.cc,
libgui/graphics/Panel.h, libgui/graphics/PopupMenuControl.cc,
libgui/graphics/PopupMenuControl.h, libgui/graphics/PushButtonControl.cc,
libgui/graphics/PushButtonControl.h, libgui/graphics/PushTool.cc,
libgui/graphics/PushTool.h, libgui/graphics/QtHandlesUtils.cc,
libgui/graphics/QtHandlesUtils.h, libgui/graphics/RadioButtonControl.cc,
libgui/graphics/RadioButtonControl.h, libgui/graphics/SliderControl.cc,
libgui/graphics/SliderControl.h, libgui/graphics/TextControl.cc,
libgui/graphics/TextControl.h, libgui/graphics/TextEdit.cc,
libgui/graphics/TextEdit.h, libgui/graphics/ToggleButtonControl.cc,
libgui/graphics/ToggleButtonControl.h, libgui/graphics/ToggleTool.cc,
libgui/graphics/ToggleTool.h, libgui/graphics/ToolBar.cc,
libgui/graphics/ToolBar.h, libgui/graphics/ToolBarButton.cc,
libgui/graphics/ToolBarButton.h, libgui/graphics/__init_qt__.cc,
libgui/graphics/__init_qt__.h, libgui/graphics/annotation-dialog.cc,
libgui/graphics/annotation-dialog.h, libgui/graphics/gl-select.cc,
libgui/graphics/gl-select.h, libgui/src/color-picker.cc,
libgui/src/color-picker.h, libgui/src/dialog.cc, libgui/src/dialog.h,
libgui/src/documentation-dock-widget.cc,
libgui/src/documentation-dock-widget.h, libgui/src/files-dock-widget.cc,
libgui/src/files-dock-widget.h, libgui/src/find-files-dialog.cc,
libgui/src/find-files-dialog.h, libgui/src/find-files-model.cc,
libgui/src/find-files-model.h, libgui/src/history-dock-widget.cc,
libgui/src/history-dock-widget.h, libgui/src/liboctgui-build-info.h,
libgui/src/liboctgui-build-info.in.cc,
libgui/src/m-editor/file-editor-interface.h,
libgui/src/m-editor/file-editor-tab.cc, libgui/src/m-editor/file-editor-tab.h,
libgui/src/m-editor/file-editor.cc, libgui/src/m-editor/file-editor.h,
libgui/src/m-editor/find-dialog.cc, libgui/src/m-editor/find-dialog.h,
libgui/src/m-editor/marker.cc, libgui/src/m-editor/marker.h,
libgui/src/m-editor/octave-qscintilla.cc,
libgui/src/m-editor/octave-qscintilla.h,
libgui/src/m-editor/octave-txt-lexer.cc,
libgui/src/m-editor/octave-txt-lexer.h, libgui/src/main-window.cc,
libgui/src/main-window.h, libgui/src/octave-cmd.cc, libgui/src/octave-cmd.h,
libgui/src/octave-dock-widget.cc, libgui/src/octave-dock-widget.h,
libgui/src/octave-gui.cc, libgui/src/octave-gui.h,
libgui/src/octave-interpreter.cc, libgui/src/octave-interpreter.h,
libgui/src/octave-qt-link.cc, libgui/src/octave-qt-link.h,
libgui/src/qtinfo/parser.cc, libgui/src/qtinfo/parser.h,
libgui/src/qtinfo/webinfo.cc, libgui/src/qtinfo/webinfo.h,
libgui/src/resource-manager.cc, libgui/src/resource-manager.h,
libgui/src/settings-dialog.cc, libgui/src/settings-dialog.h,
libgui/src/shortcut-manager.cc, libgui/src/shortcut-manager.h,
libgui/src/terminal-dock-widget.cc, libgui/src/terminal-dock-widget.h,
libgui/src/thread-manager.cc, libgui/src/thread-manager.h,
libgui/src/welcome-wizard.cc, libgui/src/welcome-wizard.h,
libgui/src/workspace-model.cc, libgui/src/workspace-model.h,
libgui/src/workspace-view.cc, libgui/src/workspace-view.h,
libinterp/build-env.h, libinterp/build-env.in.cc, libinterp/builtins.h,
libinterp/corefcn/Cell.cc, libinterp/corefcn/Cell.h,
libinterp/corefcn/__contourc__.cc, libinterp/corefcn/__dsearchn__.cc,
libinterp/corefcn/__ichol__.cc, libinterp/corefcn/__ilu__.cc,
libinterp/corefcn/__lin_interpn__.cc, libinterp/corefcn/__luinc__.cc,
libinterp/corefcn/__magick_read__.cc, libinterp/corefcn/__pchip_deriv__.cc,
libinterp/corefcn/__qp__.cc, libinterp/corefcn/balance.cc,
libinterp/corefcn/base-text-renderer.h, libinterp/corefcn/besselj.cc,
libinterp/corefcn/betainc.cc, libinterp/corefcn/bitfcns.cc,
libinterp/corefcn/bsxfun.cc, libinterp/corefcn/c-file-ptr-stream.cc,
libinterp/corefcn/c-file-ptr-stream.h, libinterp/corefcn/call-stack.cc,
libinterp/corefcn/call-stack.h, libinterp/corefcn/cdisplay.c,
libinterp/corefcn/cdisplay.h, libinterp/corefcn/cellfun.cc,
libinterp/corefcn/coct-hdf5-types.c, libinterp/corefcn/colloc.cc,
libinterp/corefcn/comment-list.cc, libinterp/corefcn/comment-list.h,
libinterp/corefcn/conv2.cc, libinterp/corefcn/daspk.cc,
libinterp/corefcn/dasrt.cc, libinterp/corefcn/dassl.cc,
libinterp/corefcn/data.cc, libinterp/corefcn/data.h,
libinterp/corefcn/debug.cc, libinterp/corefcn/debug.h,
libinterp/corefcn/defaults.cc, libinterp/corefcn/defaults.in.h,
libinterp/corefcn/defun-dld.h, libinterp/corefcn/defun-int.h,
libinterp/corefcn/defun.cc, libinterp/corefcn/defun.h,
libinterp/corefcn/det.cc, libinterp/corefcn/dirfns.cc,
libinterp/corefcn/dirfns.h, libinterp/corefcn/display.cc,
libinterp/corefcn/display.h, libinterp/corefcn/dlmread.cc,
libinterp/corefcn/dot.cc, libinterp/corefcn/dynamic-ld.cc,
libinterp/corefcn/dynamic-ld.h, libinterp/corefcn/eig.cc,
libinterp/corefcn/ellipj.cc, libinterp/corefcn/error.cc,
libinterp/corefcn/error.h, libinterp/corefcn/errwarn.cc,
libinterp/corefcn/errwarn.h, libinterp/corefcn/event-queue.cc,
libinterp/corefcn/event-queue.h, libinterp/corefcn/fft.cc,
libinterp/corefcn/fft2.cc, libinterp/corefcn/fftn.cc,
libinterp/corefcn/file-io.cc, libinterp/corefcn/file-io.h,
libinterp/corefcn/filter.cc, libinterp/corefcn/find.cc,
libinterp/corefcn/ft-text-renderer.cc, libinterp/corefcn/ft-text-renderer.h,
libinterp/corefcn/gammainc.cc, libinterp/corefcn/gcd.cc,
libinterp/corefcn/getgrent.cc, libinterp/corefcn/getpwent.cc,
libinterp/corefcn/getrusage.cc, libinterp/corefcn/givens.cc,
libinterp/corefcn/gl-render.cc, libinterp/corefcn/gl-render.h,
libinterp/corefcn/gl2ps-print.cc, libinterp/corefcn/gl2ps-print.h,
libinterp/corefcn/graphics.cc, libinterp/corefcn/graphics.in.h,
libinterp/corefcn/gripes.cc, libinterp/corefcn/gripes.h,
libinterp/corefcn/gsvd.cc, libinterp/corefcn/hash.cc,
libinterp/corefcn/help.cc, libinterp/corefcn/help.h, libinterp/corefcn/hess.cc,
libinterp/corefcn/hex2num.cc, libinterp/corefcn/hook-fcn.cc,
libinterp/corefcn/hook-fcn.h, libinterp/corefcn/input.cc,
libinterp/corefcn/input.h, libinterp/corefcn/interpreter.cc,
libinterp/corefcn/interpreter.h, libinterp/corefcn/inv.cc,
libinterp/corefcn/jit-ir.cc, libinterp/corefcn/jit-ir.h,
libinterp/corefcn/jit-typeinfo.cc, libinterp/corefcn/jit-typeinfo.h,
libinterp/corefcn/jit-util.cc, libinterp/corefcn/jit-util.h,
libinterp/corefcn/kron.cc, libinterp/corefcn/load-path.cc,
libinterp/corefcn/load-path.h, libinterp/corefcn/load-save.cc,
libinterp/corefcn/load-save.h, libinterp/corefcn/lookup.cc,
libinterp/corefcn/ls-ascii-helper.cc, libinterp/corefcn/ls-ascii-helper.h,
libinterp/corefcn/ls-hdf5.cc, libinterp/corefcn/ls-hdf5.h,
libinterp/corefcn/ls-mat-ascii.cc, libinterp/corefcn/ls-mat-ascii.h,
libinterp/corefcn/ls-mat4.cc, libinterp/corefcn/ls-mat4.h,
libinterp/corefcn/ls-mat5.cc, libinterp/corefcn/ls-mat5.h,
libinterp/corefcn/ls-oct-binary.cc, libinterp/corefcn/ls-oct-binary.h,
libinterp/corefcn/ls-oct-text.cc, libinterp/corefcn/ls-oct-text.h,
libinterp/corefcn/ls-utils.cc, libinterp/corefcn/ls-utils.h,
libinterp/corefcn/lsode.cc, libinterp/corefcn/lu.cc,
libinterp/corefcn/mappers.cc, libinterp/corefcn/matrix_type.cc,
libinterp/corefcn/max.cc, libinterp/corefcn/mex.cc, libinterp/corefcn/mex.h,
libinterp/corefcn/mexproto.h, libinterp/corefcn/mgorth.cc,
libinterp/corefcn/mxarray.in.h, libinterp/corefcn/nproc.cc,
libinterp/corefcn/oct-errno.h, libinterp/corefcn/oct-errno.in.cc,
libinterp/corefcn/oct-fstrm.cc, libinterp/corefcn/oct-fstrm.h,
libinterp/corefcn/oct-handle.h, libinterp/corefcn/oct-hdf5-types.cc,
libinterp/corefcn/oct-hdf5-types.h, libinterp/corefcn/oct-hdf5.h,
libinterp/corefcn/oct-hist.cc, libinterp/corefcn/oct-hist.h,
libinterp/corefcn/oct-iostrm.cc, libinterp/corefcn/oct-iostrm.h,
libinterp/corefcn/oct-lvalue.cc, libinterp/corefcn/oct-lvalue.h,
libinterp/corefcn/oct-map.cc, libinterp/corefcn/oct-map.h,
libinterp/corefcn/oct-obj.h, libinterp/corefcn/oct-opengl.h,
libinterp/corefcn/oct-prcstrm.cc, libinterp/corefcn/oct-prcstrm.h,
libinterp/corefcn/oct-procbuf.cc, libinterp/corefcn/oct-procbuf.h,
libinterp/corefcn/oct-stdstrm.h, libinterp/corefcn/oct-stream.cc,
libinterp/corefcn/oct-stream.h, libinterp/corefcn/oct-strstrm.cc,
libinterp/corefcn/oct-strstrm.h, libinterp/corefcn/oct-tex-lexer.in.ll,
libinterp/corefcn/oct-tex-parser.in.yy, libinterp/corefcn/oct.h,
libinterp/corefcn/octave-default-image.h, libinterp/corefcn/octave-link.cc,
libinterp/corefcn/octave-link.h,
libinterp/corefcn/octave-preserve-stream-state.h,
libinterp/corefcn/ordschur.cc, libinterp/corefcn/pager.cc,
libinterp/corefcn/pager.h, libinterp/corefcn/pinv.cc,
libinterp/corefcn/pr-output.cc, libinterp/corefcn/pr-output.h,
libinterp/corefcn/procstream.cc, libinterp/corefcn/procstream.h,
libinterp/corefcn/profiler.cc, libinterp/corefcn/profiler.h,
libinterp/corefcn/psi.cc, libinterp/corefcn/pt-jit.cc,
libinterp/corefcn/pt-jit.h, libinterp/corefcn/quad.cc,
libinterp/corefcn/quadcc.cc, libinterp/corefcn/qz.cc,
libinterp/corefcn/rand.cc, libinterp/corefcn/rcond.cc,
libinterp/corefcn/regexp.cc, libinterp/corefcn/schur.cc,
libinterp/corefcn/sighandlers.cc, libinterp/corefcn/sighandlers.h,
libinterp/corefcn/sparse-xdiv.cc, libinterp/corefcn/sparse-xdiv.h,
libinterp/corefcn/sparse-xpow.cc, libinterp/corefcn/sparse-xpow.h,
libinterp/corefcn/sparse.cc, libinterp/corefcn/spparms.cc,
libinterp/corefcn/sqrtm.cc, libinterp/corefcn/str2double.cc,
libinterp/corefcn/strfind.cc, libinterp/corefcn/strfns.cc,
libinterp/corefcn/sub2ind.cc, libinterp/corefcn/svd.cc,
libinterp/corefcn/sylvester.cc, libinterp/corefcn/symtab.cc,
libinterp/corefcn/symtab.h, libinterp/corefcn/syscalls.cc,
libinterp/corefcn/sysdep.cc, libinterp/corefcn/sysdep.h,
libinterp/corefcn/text-renderer.cc, libinterp/corefcn/text-renderer.h,
libinterp/corefcn/time.cc, libinterp/corefcn/toplev.cc,
libinterp/corefcn/toplev.h, libinterp/corefcn/tril.cc,
libinterp/corefcn/tsearch.cc, libinterp/corefcn/txt-eng.cc,
libinterp/corefcn/txt-eng.h, libinterp/corefcn/typecast.cc,
libinterp/corefcn/urlwrite.cc, libinterp/corefcn/utils.cc,
libinterp/corefcn/utils.h, libinterp/corefcn/variables.cc,
libinterp/corefcn/variables.h, libinterp/corefcn/workspace-element.h,
libinterp/corefcn/xdiv.cc, libinterp/corefcn/xdiv.h,
libinterp/corefcn/xnorm.cc, libinterp/corefcn/xnorm.h,
libinterp/corefcn/xpow.cc, libinterp/corefcn/xpow.h,
libinterp/corefcn/zfstream.cc, libinterp/corefcn/zfstream.h,
libinterp/deprecated-config.h, libinterp/dldfcn/__delaunayn__.cc,
libinterp/dldfcn/__eigs__.cc, libinterp/dldfcn/__fltk_uigetfile__.cc,
libinterp/dldfcn/__glpk__.cc, libinterp/dldfcn/__init_fltk__.cc,
libinterp/dldfcn/__init_gnuplot__.cc, libinterp/dldfcn/__osmesa_print__.cc,
libinterp/dldfcn/__voronoi__.cc, libinterp/dldfcn/amd.cc,
libinterp/dldfcn/audiodevinfo.cc, libinterp/dldfcn/audioread.cc,
libinterp/dldfcn/ccolamd.cc, libinterp/dldfcn/chol.cc,
libinterp/dldfcn/colamd.cc, libinterp/dldfcn/convhulln.cc,
libinterp/dldfcn/dmperm.cc, libinterp/dldfcn/fftw.cc, libinterp/dldfcn/gzip.cc,
libinterp/dldfcn/oct-qhull.h, libinterp/dldfcn/qr.cc,
libinterp/dldfcn/symbfact.cc, libinterp/dldfcn/symrcm.cc, libinterp/gendoc.pl,
libinterp/genprops.awk, libinterp/liboctinterp-build-info.h,
libinterp/liboctinterp-build-info.in.cc, libinterp/mk-errno-list,
libinterp/mk-pkg-add, libinterp/mkops, libinterp/octave-value/ov-base-diag.cc,
libinterp/octave-value/ov-base-diag.h, libinterp/octave-value/ov-base-int.cc,
libinterp/octave-value/ov-base-int.h, libinterp/octave-value/ov-base-mat.cc,
libinterp/octave-value/ov-base-mat.h, libinterp/octave-value/ov-base-scalar.cc,
libinterp/octave-value/ov-base-scalar.h,
libinterp/octave-value/ov-base-sparse.cc,
libinterp/octave-value/ov-base-sparse.h, libinterp/octave-value/ov-base.cc,
libinterp/octave-value/ov-base.h, libinterp/octave-value/ov-bool-mat.cc,
libinterp/octave-value/ov-bool-mat.h, libinterp/octave-value/ov-bool-sparse.cc,
libinterp/octave-value/ov-bool-sparse.h, libinterp/octave-value/ov-bool.cc,
libinterp/octave-value/ov-bool.h, libinterp/octave-value/ov-builtin.cc,
libinterp/octave-value/ov-builtin.h, libinterp/octave-value/ov-cell.cc,
libinterp/octave-value/ov-cell.h, libinterp/octave-value/ov-ch-mat.cc,
libinterp/octave-value/ov-ch-mat.h, libinterp/octave-value/ov-class.cc,
libinterp/octave-value/ov-class.h, libinterp/octave-value/ov-classdef.cc,
libinterp/octave-value/ov-classdef.h, libinterp/octave-value/ov-colon.cc,
libinterp/octave-value/ov-colon.h, libinterp/octave-value/ov-complex.cc,
libinterp/octave-value/ov-complex.h, libinterp/octave-value/ov-cs-list.cc,
libinterp/octave-value/ov-cs-list.h, libinterp/octave-value/ov-cx-diag.cc,
libinterp/octave-value/ov-cx-diag.h, libinterp/octave-value/ov-cx-mat.cc,
libinterp/octave-value/ov-cx-mat.h, libinterp/octave-value/ov-cx-sparse.cc,
libinterp/octave-value/ov-cx-sparse.h, libinterp/octave-value/ov-dld-fcn.cc,
libinterp/octave-value/ov-dld-fcn.h, libinterp/octave-value/ov-fcn-handle.cc,
libinterp/octave-value/ov-fcn-handle.h,
libinterp/octave-value/ov-fcn-inline.cc,
libinterp/octave-value/ov-fcn-inline.h, libinterp/octave-value/ov-fcn.cc,
libinterp/octave-value/ov-fcn.h, libinterp/octave-value/ov-float.cc,
libinterp/octave-value/ov-float.h, libinterp/octave-value/ov-flt-complex.cc,
libinterp/octave-value/ov-flt-complex.h,
libinterp/octave-value/ov-flt-cx-diag.cc,
libinterp/octave-value/ov-flt-cx-diag.h,
libinterp/octave-value/ov-flt-cx-mat.cc,
libinterp/octave-value/ov-flt-cx-mat.h,
libinterp/octave-value/ov-flt-re-diag.cc,
libinterp/octave-value/ov-flt-re-diag.h,
libinterp/octave-value/ov-flt-re-mat.cc,
libinterp/octave-value/ov-flt-re-mat.h, libinterp/octave-value/ov-int-traits.h,
libinterp/octave-value/ov-int16.cc, libinterp/octave-value/ov-int16.h,
libinterp/octave-value/ov-int32.cc, libinterp/octave-value/ov-int32.h,
libinterp/octave-value/ov-int64.cc, libinterp/octave-value/ov-int64.h,
libinterp/octave-value/ov-int8.cc, libinterp/octave-value/ov-int8.h,
libinterp/octave-value/ov-intx.h, libinterp/octave-value/ov-java.cc,
libinterp/octave-value/ov-java.h, libinterp/octave-value/ov-lazy-idx.cc,
libinterp/octave-value/ov-lazy-idx.h, libinterp/octave-value/ov-mex-fcn.cc,
libinterp/octave-value/ov-mex-fcn.h, libinterp/octave-value/ov-null-mat.cc,
libinterp/octave-value/ov-null-mat.h, libinterp/octave-value/ov-oncleanup.cc,
libinterp/octave-value/ov-oncleanup.h, libinterp/octave-value/ov-perm.cc,
libinterp/octave-value/ov-perm.h, libinterp/octave-value/ov-range.cc,
libinterp/octave-value/ov-range.h, libinterp/octave-value/ov-re-diag.cc,
libinterp/octave-value/ov-re-diag.h, libinterp/octave-value/ov-re-mat.cc,
libinterp/octave-value/ov-re-mat.h, libinterp/octave-value/ov-re-sparse.cc,
libinterp/octave-value/ov-re-sparse.h, libinterp/octave-value/ov-scalar.cc,
libinterp/octave-value/ov-scalar.h, libinterp/octave-value/ov-str-mat.cc,
libinterp/octave-value/ov-str-mat.h, libinterp/octave-value/ov-struct.cc,
libinterp/octave-value/ov-struct.h, libinterp/octave-value/ov-typeinfo.cc,
libinterp/octave-value/ov-typeinfo.h, libinterp/octave-value/ov-uint16.cc,
libinterp/octave-value/ov-uint16.h, libinterp/octave-value/ov-uint32.cc,
libinterp/octave-value/ov-uint32.h, libinterp/octave-value/ov-uint64.cc,
libinterp/octave-value/ov-uint64.h, libinterp/octave-value/ov-uint8.cc,
libinterp/octave-value/ov-uint8.h, libinterp/octave-value/ov-usr-fcn.cc,
libinterp/octave-value/ov-usr-fcn.h, libinterp/octave-value/ov.cc,
libinterp/octave-value/ov.h, libinterp/octave-value/ovl.cc,
libinterp/octave-value/ovl.h, libinterp/octave.cc, libinterp/octave.h,
libinterp/op-kw-docs, libinterp/operators/op-b-b.cc,
libinterp/operators/op-b-bm.cc, libinterp/operators/op-b-sbm.cc,
libinterp/operators/op-bm-b.cc, libinterp/operators/op-bm-bm.cc,
libinterp/operators/op-bm-sbm.cc, libinterp/operators/op-cdm-cdm.cc,
libinterp/operators/op-cdm-cm.cc, libinterp/operators/op-cdm-cs.cc,
libinterp/operators/op-cdm-dm.cc, libinterp/operators/op-cdm-m.cc,
libinterp/operators/op-cdm-s.cc, libinterp/operators/op-cell.cc,
libinterp/operators/op-chm.cc, libinterp/operators/op-class.cc,
libinterp/operators/op-cm-cdm.cc, libinterp/operators/op-cm-cm.cc,
libinterp/operators/op-cm-cs.cc, libinterp/operators/op-cm-dm.cc,
libinterp/operators/op-cm-m.cc, libinterp/operators/op-cm-pm.cc,
libinterp/operators/op-cm-s.cc, libinterp/operators/op-cm-scm.cc,
libinterp/operators/op-cm-sm.cc, libinterp/operators/op-cs-cm.cc,
libinterp/operators/op-cs-cs.cc, libinterp/operators/op-cs-m.cc,
libinterp/operators/op-cs-s.cc, libinterp/operators/op-cs-scm.cc,
libinterp/operators/op-cs-sm.cc, libinterp/operators/op-dm-cdm.cc,
libinterp/operators/op-dm-cm.cc, libinterp/operators/op-dm-cs.cc,
libinterp/operators/op-dm-dm.cc, libinterp/operators/op-dm-m.cc,
libinterp/operators/op-dm-s.cc, libinterp/operators/op-dm-scm.cc,
libinterp/operators/op-dm-sm.cc, libinterp/operators/op-dm-template.cc,
libinterp/operators/op-dms-template.cc, libinterp/operators/op-fcdm-fcdm.cc,
libinterp/operators/op-fcdm-fcm.cc, libinterp/operators/op-fcdm-fcs.cc,
libinterp/operators/op-fcdm-fdm.cc, libinterp/operators/op-fcdm-fm.cc,
libinterp/operators/op-fcdm-fs.cc, libinterp/operators/op-fcm-fcdm.cc,
libinterp/operators/op-fcm-fcm.cc, libinterp/operators/op-fcm-fcs.cc,
libinterp/operators/op-fcm-fdm.cc, libinterp/operators/op-fcm-fm.cc,
libinterp/operators/op-fcm-fs.cc, libinterp/operators/op-fcm-pm.cc,
libinterp/operators/op-fcn.cc, libinterp/operators/op-fcs-fcm.cc,
libinterp/operators/op-fcs-fcs.cc, libinterp/operators/op-fcs-fm.cc,
libinterp/operators/op-fcs-fs.cc, libinterp/operators/op-fdm-fcdm.cc,
libinterp/operators/op-fdm-fcm.cc, libinterp/operators/op-fdm-fcs.cc,
libinterp/operators/op-fdm-fdm.cc, libinterp/operators/op-fdm-fm.cc,
libinterp/operators/op-fdm-fs.cc, libinterp/operators/op-fm-fcdm.cc,
libinterp/operators/op-fm-fcm.cc, libinterp/operators/op-fm-fcs.cc,
libinterp/operators/op-fm-fdm.cc, libinterp/operators/op-fm-fm.cc,
libinterp/operators/op-fm-fs.cc, libinterp/operators/op-fm-pm.cc,
libinterp/operators/op-fs-fcm.cc, libinterp/operators/op-fs-fcs.cc,
libinterp/operators/op-fs-fm.cc, libinterp/operators/op-fs-fs.cc,
libinterp/operators/op-i16-i16.cc, libinterp/operators/op-i32-i32.cc,
libinterp/operators/op-i64-i64.cc, libinterp/operators/op-i8-i8.cc,
libinterp/operators/op-int-concat.cc, libinterp/operators/op-int.h,
libinterp/operators/op-m-cdm.cc, libinterp/operators/op-m-cm.cc,
libinterp/operators/op-m-cs.cc, libinterp/operators/op-m-dm.cc,
libinterp/operators/op-m-m.cc, libinterp/operators/op-m-pm.cc,
libinterp/operators/op-m-s.cc, libinterp/operators/op-m-scm.cc,
libinterp/operators/op-m-sm.cc, libinterp/operators/op-pm-cm.cc,
libinterp/operators/op-pm-fcm.cc, libinterp/operators/op-pm-fm.cc,
libinterp/operators/op-pm-m.cc, libinterp/operators/op-pm-pm.cc,
libinterp/operators/op-pm-scm.cc, libinterp/operators/op-pm-sm.cc,
libinterp/operators/op-pm-template.cc, libinterp/operators/op-range.cc,
libinterp/operators/op-s-cm.cc, libinterp/operators/op-s-cs.cc,
libinterp/operators/op-s-m.cc, libinterp/operators/op-s-s.cc,
libinterp/operators/op-s-scm.cc, libinterp/operators/op-s-sm.cc,
libinterp/operators/op-sbm-b.cc, libinterp/operators/op-sbm-bm.cc,
libinterp/operators/op-sbm-sbm.cc, libinterp/operators/op-scm-cm.cc,
libinterp/operators/op-scm-cs.cc, libinterp/operators/op-scm-m.cc,
libinterp/operators/op-scm-s.cc, libinterp/operators/op-scm-scm.cc,
libinterp/operators/op-scm-sm.cc, libinterp/operators/op-sm-cm.cc,
libinterp/operators/op-sm-cs.cc, libinterp/operators/op-sm-m.cc,
libinterp/operators/op-sm-s.cc, libinterp/operators/op-sm-scm.cc,
libinterp/operators/op-sm-sm.cc, libinterp/operators/op-str-m.cc,
libinterp/operators/op-str-s.cc, libinterp/operators/op-str-str.cc,
libinterp/operators/op-struct.cc, libinterp/operators/op-ui16-ui16.cc,
libinterp/operators/op-ui32-ui32.cc, libinterp/operators/op-ui64-ui64.cc,
libinterp/operators/op-ui8-ui8.cc, libinterp/operators/ops.h,
libinterp/options-usage.h, libinterp/parse-tree/lex.h,
libinterp/parse-tree/lex.ll, libinterp/parse-tree/oct-parse.in.yy,
libinterp/parse-tree/octave.gperf, libinterp/parse-tree/parse.h,
libinterp/parse-tree/pt-all.h, libinterp/parse-tree/pt-arg-list.cc,
libinterp/parse-tree/pt-arg-list.h, libinterp/parse-tree/pt-array-list.cc,
libinterp/parse-tree/pt-array-list.h, libinterp/parse-tree/pt-assign.cc,
libinterp/parse-tree/pt-assign.h, libinterp/parse-tree/pt-binop.cc,
libinterp/parse-tree/pt-binop.h, libinterp/parse-tree/pt-bp.cc,
libinterp/parse-tree/pt-bp.h, libinterp/parse-tree/pt-cbinop.cc,
libinterp/parse-tree/pt-cbinop.h, libinterp/parse-tree/pt-cell.cc,
libinterp/parse-tree/pt-cell.h, libinterp/parse-tree/pt-check.cc,
libinterp/parse-tree/pt-check.h, libinterp/parse-tree/pt-classdef.cc,
libinterp/parse-tree/pt-classdef.h, libinterp/parse-tree/pt-cmd.cc,
libinterp/parse-tree/pt-cmd.h, libinterp/parse-tree/pt-colon.cc,
libinterp/parse-tree/pt-colon.h, libinterp/parse-tree/pt-const.cc,
libinterp/parse-tree/pt-const.h, libinterp/parse-tree/pt-decl.cc,
libinterp/parse-tree/pt-decl.h, libinterp/parse-tree/pt-eval.cc,
libinterp/parse-tree/pt-eval.h, libinterp/parse-tree/pt-except.cc,
libinterp/parse-tree/pt-except.h, libinterp/parse-tree/pt-exp.cc,
libinterp/parse-tree/pt-exp.h, libinterp/parse-tree/pt-fcn-handle.cc,
libinterp/parse-tree/pt-fcn-handle.h, libinterp/parse-tree/pt-funcall.cc,
libinterp/parse-tree/pt-funcall.h, libinterp/parse-tree/pt-id.cc,
libinterp/parse-tree/pt-id.h, libinterp/parse-tree/pt-idx.cc,
libinterp/parse-tree/pt-idx.h, libinterp/parse-tree/pt-jump.cc,
libinterp/parse-tree/pt-jump.h, libinterp/parse-tree/pt-loop.cc,
libinterp/parse-tree/pt-loop.h, libinterp/parse-tree/pt-mat.cc,
libinterp/parse-tree/pt-mat.h, libinterp/parse-tree/pt-misc.cc,
libinterp/parse-tree/pt-misc.h, libinterp/parse-tree/pt-pr-code.cc,
libinterp/parse-tree/pt-pr-code.h, libinterp/parse-tree/pt-select.cc,
libinterp/parse-tree/pt-select.h, libinterp/parse-tree/pt-stmt.cc,
libinterp/parse-tree/pt-stmt.h, libinterp/parse-tree/pt-unop.cc,
libinterp/parse-tree/pt-unop.h, libinterp/parse-tree/pt-walk.h,
libinterp/parse-tree/pt.cc, libinterp/parse-tree/pt.h,
libinterp/parse-tree/token.cc, libinterp/parse-tree/token.h,
libinterp/template-inst/Array-jit.cc, libinterp/template-inst/Array-tc.cc,
libinterp/version.cc, libinterp/version.in.h, liboctave/array/Array-C.cc,
liboctave/array/Array-b.cc, liboctave/array/Array-ch.cc,
liboctave/array/Array-d.cc, liboctave/array/Array-f.cc,
liboctave/array/Array-fC.cc, liboctave/array/Array-i.cc,
liboctave/array/Array-idx-vec.cc, liboctave/array/Array-s.cc,
liboctave/array/Array-str.cc, liboctave/array/Array-util.cc,
liboctave/array/Array-util.h, liboctave/array/Array-voidp.cc,
liboctave/array/Array.cc, liboctave/array/Array.h,
liboctave/array/CColVector.cc, liboctave/array/CColVector.h,
liboctave/array/CDiagMatrix.cc, liboctave/array/CDiagMatrix.h,
liboctave/array/CMatrix.cc, liboctave/array/CMatrix.h,
liboctave/array/CNDArray.cc, liboctave/array/CNDArray.h,
liboctave/array/CRowVector.cc, liboctave/array/CRowVector.h,
liboctave/array/CSparse.cc, liboctave/array/CSparse.h,
liboctave/array/DiagArray2.cc, liboctave/array/DiagArray2.h,
liboctave/array/MArray-C.cc, liboctave/array/MArray-d.cc,
liboctave/array/MArray-f.cc, liboctave/array/MArray-fC.cc,
liboctave/array/MArray-i.cc, liboctave/array/MArray-s.cc,
liboctave/array/MArray.cc, liboctave/array/MArray.h,
liboctave/array/MDiagArray2.cc, liboctave/array/MDiagArray2.h,
liboctave/array/MSparse-C.cc, liboctave/array/MSparse-d.cc,
liboctave/array/MSparse.cc, liboctave/array/MSparse.h,
liboctave/array/Matrix.h, liboctave/array/MatrixType.cc,
liboctave/array/MatrixType.h, liboctave/array/PermMatrix.cc,
liboctave/array/PermMatrix.h, liboctave/array/Range.cc,
liboctave/array/Range.h, liboctave/array/Sparse-C.cc,
liboctave/array/Sparse-b.cc, liboctave/array/Sparse-d.cc,
liboctave/array/Sparse.cc, liboctave/array/Sparse.h,
liboctave/array/boolMatrix.cc, liboctave/array/boolMatrix.h,
liboctave/array/boolNDArray.cc, liboctave/array/boolNDArray.h,
liboctave/array/boolSparse.cc, liboctave/array/boolSparse.h,
liboctave/array/chMatrix.cc, liboctave/array/chMatrix.h,
liboctave/array/chNDArray.cc, liboctave/array/chNDArray.h,
liboctave/array/dColVector.cc, liboctave/array/dColVector.h,
liboctave/array/dDiagMatrix.cc, liboctave/array/dDiagMatrix.h,
liboctave/array/dMatrix.cc, liboctave/array/dMatrix.h,
liboctave/array/dNDArray.cc, liboctave/array/dNDArray.h,
liboctave/array/dRowVector.cc, liboctave/array/dRowVector.h,
liboctave/array/dSparse.cc, liboctave/array/dSparse.h,
liboctave/array/dim-vector.cc, liboctave/array/dim-vector.h,
liboctave/array/fCColVector.cc, liboctave/array/fCColVector.h,
liboctave/array/fCDiagMatrix.cc, liboctave/array/fCDiagMatrix.h,
liboctave/array/fCMatrix.cc, liboctave/array/fCMatrix.h,
liboctave/array/fCNDArray.cc, liboctave/array/fCNDArray.h,
liboctave/array/fCRowVector.cc, liboctave/array/fCRowVector.h,
liboctave/array/fColVector.cc, liboctave/array/fColVector.h,
liboctave/array/fDiagMatrix.cc, liboctave/array/fDiagMatrix.h,
liboctave/array/fMatrix.cc, liboctave/array/fMatrix.h,
liboctave/array/fNDArray.cc, liboctave/array/fNDArray.h,
liboctave/array/fRowVector.cc, liboctave/array/fRowVector.h,
liboctave/array/idx-vector.cc, liboctave/array/idx-vector.h,
liboctave/array/int16NDArray.cc, liboctave/array/int16NDArray.h,
liboctave/array/int32NDArray.cc, liboctave/array/int32NDArray.h,
liboctave/array/int64NDArray.cc, liboctave/array/int64NDArray.h,
liboctave/array/int8NDArray.cc, liboctave/array/int8NDArray.h,
liboctave/array/intNDArray.cc, liboctave/array/intNDArray.h,
liboctave/array/uint16NDArray.cc, liboctave/array/uint16NDArray.h,
liboctave/array/uint32NDArray.cc, liboctave/array/uint32NDArray.h,
liboctave/array/uint64NDArray.cc, liboctave/array/uint64NDArray.h,
liboctave/array/uint8NDArray.cc, liboctave/array/uint8NDArray.h,
liboctave/cruft/misc/blaswrap.c, liboctave/cruft/misc/cquit.c,
liboctave/cruft/misc/f77-extern.cc, liboctave/cruft/misc/f77-fcn.c,
liboctave/cruft/misc/f77-fcn.h, liboctave/cruft/misc/lo-error.c,
liboctave/cruft/misc/lo-error.h, liboctave/cruft/misc/quit.cc,
liboctave/cruft/misc/quit.h, liboctave/liboctave-build-info.h,
liboctave/liboctave-build-info.in.cc, liboctave/numeric/CollocWt.cc,
liboctave/numeric/CollocWt.h, liboctave/numeric/DAE.h,
liboctave/numeric/DAEFunc.h, liboctave/numeric/DAERT.h,
liboctave/numeric/DAERTFunc.h, liboctave/numeric/DASPK-opts.in,
liboctave/numeric/DASPK.cc, liboctave/numeric/DASPK.h,
liboctave/numeric/DASRT-opts.in, liboctave/numeric/DASRT.cc,
liboctave/numeric/DASRT.h, liboctave/numeric/DASSL-opts.in,
liboctave/numeric/DASSL.cc, liboctave/numeric/DASSL.h, liboctave/numeric/DET.h,
liboctave/numeric/EIG.cc, liboctave/numeric/EIG.h,
liboctave/numeric/LSODE-opts.in, liboctave/numeric/LSODE.cc,
liboctave/numeric/LSODE.h, liboctave/numeric/ODE.h,
liboctave/numeric/ODEFunc.h, liboctave/numeric/ODES.cc,
liboctave/numeric/ODES.h, liboctave/numeric/ODESFunc.h,
liboctave/numeric/Quad-opts.in, liboctave/numeric/Quad.cc,
liboctave/numeric/Quad.h, liboctave/numeric/aepbalance.cc,
liboctave/numeric/aepbalance.h, liboctave/numeric/base-dae.h,
liboctave/numeric/base-de.h, liboctave/numeric/base-min.h,
liboctave/numeric/bsxfun-decl.h, liboctave/numeric/bsxfun-defs.cc,
liboctave/numeric/bsxfun.h, liboctave/numeric/chol.cc,
liboctave/numeric/chol.h, liboctave/numeric/eigs-base.cc,
liboctave/numeric/eigs-base.h, liboctave/numeric/fEIG.cc,
liboctave/numeric/fEIG.h, liboctave/numeric/gepbalance.cc,
liboctave/numeric/gepbalance.h, liboctave/numeric/gsvd.cc,
liboctave/numeric/gsvd.h, liboctave/numeric/hess.cc, liboctave/numeric/hess.h,
liboctave/numeric/lo-amos-proto.h, liboctave/numeric/lo-arpack-proto.h,
liboctave/numeric/lo-blas-proto.h, liboctave/numeric/lo-fftpack-proto.h,
liboctave/numeric/lo-lapack-proto.h, liboctave/numeric/lo-mappers.cc,
liboctave/numeric/lo-mappers.h, liboctave/numeric/lo-qrupdate-proto.h,
liboctave/numeric/lo-ranlib-proto.h, liboctave/numeric/lo-slatec-proto.h,
liboctave/numeric/lo-specfun.cc, liboctave/numeric/lo-specfun.h,
liboctave/numeric/lu.cc, liboctave/numeric/lu.h,
liboctave/numeric/oct-convn.cc, liboctave/numeric/oct-convn.h,
liboctave/numeric/oct-fftw.cc, liboctave/numeric/oct-fftw.h,
liboctave/numeric/oct-norm.cc, liboctave/numeric/oct-norm.h,
liboctave/numeric/oct-rand.cc, liboctave/numeric/oct-rand.h,
liboctave/numeric/oct-spparms.cc, liboctave/numeric/oct-spparms.h,
liboctave/numeric/qr.cc, liboctave/numeric/qr.h, liboctave/numeric/qrp.cc,
liboctave/numeric/qrp.h, liboctave/numeric/randgamma.cc,
liboctave/numeric/randgamma.h, liboctave/numeric/randmtzig.cc,
liboctave/numeric/randmtzig.h, liboctave/numeric/randpoisson.cc,
liboctave/numeric/randpoisson.h, liboctave/numeric/schur.cc,
liboctave/numeric/schur.h, liboctave/numeric/sparse-chol.cc,
liboctave/numeric/sparse-chol.h, liboctave/numeric/sparse-dmsolve.cc,
liboctave/numeric/sparse-dmsolve.h, liboctave/numeric/sparse-lu.cc,
liboctave/numeric/sparse-lu.h, liboctave/numeric/sparse-qr.cc,
liboctave/numeric/sparse-qr.h, liboctave/numeric/svd.cc,
liboctave/numeric/svd.h, liboctave/operators/Sparse-diag-op-defs.h,
liboctave/operators/Sparse-op-decls.h, liboctave/operators/Sparse-op-defs.h,
liboctave/operators/Sparse-perm-op-defs.h, liboctave/operators/mk-ops.awk,
liboctave/operators/mx-base.h, liboctave/operators/mx-defs.h,
liboctave/operators/mx-ext.h, liboctave/operators/mx-inlines.cc,
liboctave/operators/mx-op-decl.h, liboctave/operators/mx-op-defs.h,
liboctave/operators/mx-ops, liboctave/operators/smx-ops,
liboctave/operators/vx-ops, liboctave/system/child-list.cc,
liboctave/system/child-list.h, liboctave/system/dir-ops.cc,
liboctave/system/dir-ops.h, liboctave/system/file-ops.cc,
liboctave/system/file-ops.h, liboctave/system/file-stat.cc,
liboctave/system/file-stat.h, liboctave/system/lo-sysdep.cc,
liboctave/system/lo-sysdep.h, liboctave/system/mach-info.cc,
liboctave/system/mach-info.h, liboctave/system/oct-env.cc,
liboctave/system/oct-env.h, liboctave/system/oct-group.cc,
liboctave/system/oct-group.h, liboctave/system/oct-passwd.cc,
liboctave/system/oct-passwd.h, liboctave/system/oct-syscalls.cc,
liboctave/system/oct-syscalls.h, liboctave/system/oct-time.cc,
liboctave/system/oct-time.h, liboctave/system/oct-uname.cc,
liboctave/system/oct-uname.h, liboctave/util/action-container.h,
liboctave/util/base-list.h, liboctave/util/byte-swap.h,
liboctave/util/caseless-str.h, liboctave/util/cmd-edit.cc,
liboctave/util/cmd-edit.h, liboctave/util/cmd-hist.cc,
liboctave/util/cmd-hist.h, liboctave/util/data-conv.cc,
liboctave/util/data-conv.h, liboctave/util/f2c-main.c,
liboctave/util/functor.h, liboctave/util/glob-match.cc,
liboctave/util/glob-match.h, liboctave/util/kpse.cc, liboctave/util/kpse.h,
liboctave/util/lo-array-errwarn.cc, liboctave/util/lo-array-errwarn.h,
liboctave/util/lo-array-gripes.cc, liboctave/util/lo-array-gripes.h,
liboctave/util/lo-cutils.c, liboctave/util/lo-cutils.h,
liboctave/util/lo-hash.cc, liboctave/util/lo-hash.h, liboctave/util/lo-ieee.cc,
liboctave/util/lo-ieee.h, liboctave/util/lo-macros.h, liboctave/util/lo-math.h,
liboctave/util/lo-regexp.cc, liboctave/util/lo-regexp.h,
liboctave/util/lo-traits.h, liboctave/util/lo-utils.cc,
liboctave/util/lo-utils.h, liboctave/util/oct-alloc.h,
liboctave/util/oct-base64.cc, liboctave/util/oct-base64.h,
liboctave/util/oct-binmap.h, liboctave/util/oct-cmplx.h,
liboctave/util/oct-glob.cc, liboctave/util/oct-glob.h,
liboctave/util/oct-inttypes-fwd.h, liboctave/util/oct-inttypes.cc,
liboctave/util/oct-inttypes.h, liboctave/util/oct-locbuf.cc,
liboctave/util/oct-locbuf.h, liboctave/util/oct-mutex.cc,
liboctave/util/oct-mutex.h, liboctave/util/oct-refcount.h,
liboctave/util/oct-rl-edit.c, liboctave/util/oct-rl-edit.h,
liboctave/util/oct-rl-hist.c, liboctave/util/oct-rl-hist.h,
liboctave/util/oct-shlib.cc, liboctave/util/oct-shlib.h,
liboctave/util/oct-sort.cc, liboctave/util/oct-sort.h,
liboctave/util/oct-sparse.h, liboctave/util/oct-string.cc,
liboctave/util/oct-string.h, liboctave/util/pathsearch.cc,
liboctave/util/pathsearch.h, liboctave/util/singleton-cleanup.cc,
liboctave/util/sparse-sort.cc, liboctave/util/sparse-sort.h,
liboctave/util/sparse-util.cc, liboctave/util/sparse-util.h,
liboctave/util/str-vec.cc, liboctave/util/str-vec.h,
liboctave/util/sun-utils.h, liboctave/util/unwind-prot.cc,
liboctave/util/unwind-prot.h, liboctave/util/url-transfer.cc,
liboctave/util/url-transfer.h, liboctave/wrappers/areadlink-wrapper.c,
liboctave/wrappers/areadlink-wrapper.h,
liboctave/wrappers/async-system-wrapper.c,
liboctave/wrappers/async-system-wrapper.h,
liboctave/wrappers/base64-wrappers.c, liboctave/wrappers/base64-wrappers.h,
liboctave/wrappers/canonicalize-file-name-wrapper.c,
liboctave/wrappers/canonicalize-file-name-wrapper.h,
liboctave/wrappers/dirent-wrappers.c, liboctave/wrappers/dirent-wrappers.h,
liboctave/wrappers/fcntl-wrappers.c, liboctave/wrappers/fcntl-wrappers.h,
liboctave/wrappers/filepos-wrappers.c, liboctave/wrappers/filepos-wrappers.h,
liboctave/wrappers/fpucw-wrappers.c, liboctave/wrappers/fpucw-wrappers.h,
liboctave/wrappers/gen-tempname-wrapper.c,
liboctave/wrappers/gen-tempname-wrapper.h, liboctave/wrappers/getopt-wrapper.c,
liboctave/wrappers/getopt-wrapper.h, liboctave/wrappers/glob-wrappers.c,
liboctave/wrappers/glob-wrappers.h, liboctave/wrappers/hash-wrappers.c,
liboctave/wrappers/hash-wrappers.h, liboctave/wrappers/math-wrappers.c,
liboctave/wrappers/math-wrappers.h, liboctave/wrappers/mkostemp-wrapper.c,
liboctave/wrappers/mkostemp-wrapper.h, liboctave/wrappers/nanosleep-wrapper.c,
liboctave/wrappers/nanosleep-wrapper.h, liboctave/wrappers/nproc-wrapper.c,
liboctave/wrappers/nproc-wrapper.h, liboctave/wrappers/octave-popen2.c,
liboctave/wrappers/octave-popen2.h, liboctave/wrappers/putenv-wrapper.c,
liboctave/wrappers/putenv-wrapper.h,
liboctave/wrappers/set-program-name-wrapper.c,
liboctave/wrappers/set-program-name-wrapper.h,
liboctave/wrappers/signal-wrappers.c, liboctave/wrappers/signal-wrappers.h,
liboctave/wrappers/stat-wrappers.c, liboctave/wrappers/stat-wrappers.h,
liboctave/wrappers/strdup-wrapper.c, liboctave/wrappers/strdup-wrapper.h,
liboctave/wrappers/strftime-wrapper.c, liboctave/wrappers/strftime-wrapper.h,
liboctave/wrappers/strmode-wrapper.c, liboctave/wrappers/strmode-wrapper.h,
liboctave/wrappers/strptime-wrapper.c, liboctave/wrappers/strptime-wrapper.h,
liboctave/wrappers/time-wrappers.c, liboctave/wrappers/time-wrappers.h,
liboctave/wrappers/tmpfile-wrapper.c, liboctave/wrappers/tmpfile-wrapper.h,
liboctave/wrappers/uname-wrapper.c, liboctave/wrappers/uname-wrapper.h,
liboctave/wrappers/unistd-wrappers.c, liboctave/wrappers/unistd-wrappers.h,
liboctave/wrappers/unsetenv-wrapper.c, liboctave/wrappers/unsetenv-wrapper.h,
liboctave/wrappers/vasprintf-wrapper.c, liboctave/wrappers/vasprintf-wrapper.h,
liboctave/wrappers/wait-for-input.c, liboctave/wrappers/wait-for-input.h,
liboctave/wrappers/wait-wrappers.c, liboctave/wrappers/wait-wrappers.h,
m4/acinclude.m4, oct-conf-post.in.h, run-octave.in, scripts/@ftp/ascii.m,
scripts/@ftp/binary.m, scripts/@ftp/cd.m, scripts/@ftp/close.m,
scripts/@ftp/delete.m, scripts/@ftp/dir.m, scripts/@ftp/display.m,
scripts/@ftp/ftp.m, scripts/@ftp/loadobj.m, scripts/@ftp/mget.m,
scripts/@ftp/mkdir.m, scripts/@ftp/mput.m, scripts/@ftp/rename.m,
scripts/@ftp/rmdir.m, scripts/@ftp/saveobj.m,
scripts/audio/@audioplayer/__get_properties__.m,
scripts/audio/@audioplayer/audioplayer.m, scripts/audio/@audioplayer/display.m,
scripts/audio/@audioplayer/get.m, scripts/audio/@audioplayer/isplaying.m,
scripts/audio/@audioplayer/pause.m, scripts/audio/@audioplayer/play.m,
scripts/audio/@audioplayer/playblocking.m, scripts/audio/@audioplayer/resume.m,
scripts/audio/@audioplayer/set.m, scripts/audio/@audioplayer/stop.m,
scripts/audio/@audioplayer/subsasgn.m, scripts/audio/@audioplayer/subsref.m,
scripts/audio/@audiorecorder/__get_properties__.m,
scripts/audio/@audiorecorder/audiorecorder.m,
scripts/audio/@audiorecorder/display.m, scripts/audio/@audiorecorder/get.m,
scripts/audio/@audiorecorder/getaudiodata.m,
scripts/audio/@audiorecorder/getplayer.m,
scripts/audio/@audiorecorder/isrecording.m,
scripts/audio/@audiorecorder/pause.m, scripts/audio/@audiorecorder/play.m,
scripts/audio/@audiorecorder/record.m,
scripts/audio/@audiorecorder/recordblocking.m,
scripts/audio/@audiorecorder/resume.m, scripts/audio/@audiorecorder/set.m,
scripts/audio/@audiorecorder/stop.m, scripts/audio/@audiorecorder/subsasgn.m,
scripts/audio/@audiorecorder/subsref.m, scripts/audio/lin2mu.m,
scripts/audio/mu2lin.m, scripts/audio/record.m, scripts/audio/sound.m,
scripts/audio/soundsc.m, scripts/deprecated/bitmax.m,
scripts/deprecated/comma.m, scripts/deprecated/isstr.m,
scripts/deprecated/mahalanobis.m, scripts/deprecated/md5sum.m,
scripts/deprecated/octave_config_info.m, scripts/deprecated/onenormest.m,
scripts/deprecated/paren.m, scripts/deprecated/semicolon.m,
scripts/deprecated/sleep.m, scripts/deprecated/usleep.m,
scripts/deprecated/wavread.m, scripts/deprecated/wavwrite.m,
scripts/elfun/acosd.m, scripts/elfun/acot.m, scripts/elfun/acotd.m,
scripts/elfun/acoth.m, scripts/elfun/acsc.m, scripts/elfun/acscd.m,
scripts/elfun/acsch.m, scripts/elfun/asec.m, scripts/elfun/asecd.m,
scripts/elfun/asech.m, scripts/elfun/asind.m, scripts/elfun/atan2d.m,
scripts/elfun/atand.m, scripts/elfun/cosd.m, scripts/elfun/cot.m,
scripts/elfun/cotd.m, scripts/elfun/coth.m, scripts/elfun/csc.m,
scripts/elfun/cscd.m, scripts/elfun/csch.m, scripts/elfun/sec.m,
scripts/elfun/secd.m, scripts/elfun/sech.m, scripts/elfun/sind.m,
scripts/elfun/tand.m, scripts/general/accumarray.m, scripts/general/accumdim.m,
scripts/general/bincoeff.m, scripts/general/bitcmp.m, scripts/general/bitget.m,
scripts/general/bitset.m, scripts/general/blkdiag.m,
scripts/general/cart2pol.m, scripts/general/cart2sph.m,
scripts/general/cell2mat.m, scripts/general/celldisp.m, scripts/general/chop.m,
scripts/general/circshift.m, scripts/general/common_size.m,
scripts/general/cplxpair.m, scripts/general/cumtrapz.m, scripts/general/curl.m,
scripts/general/dblquad.m, scripts/general/deal.m, scripts/general/deg2rad.m,
scripts/general/del2.m, scripts/general/display.m,
scripts/general/divergence.m, scripts/general/fieldnames.m,
scripts/general/flip.m, scripts/general/flipdim.m, scripts/general/fliplr.m,
scripts/general/flipud.m, scripts/general/grabcode.m,
scripts/general/gradient.m, scripts/general/idivide.m,
scripts/general/inputParser.m, scripts/general/int2str.m,
scripts/general/interp1.m, scripts/general/interp2.m,
scripts/general/interp3.m, scripts/general/interpft.m,
scripts/general/interpn.m, scripts/general/isdir.m, scripts/general/isequal.m,
scripts/general/isequaln.m, scripts/general/loadobj.m,
scripts/general/logspace.m, scripts/general/methods.m,
scripts/general/nargchk.m, scripts/general/narginchk.m,
scripts/general/nargoutchk.m, scripts/general/nextpow2.m,
scripts/general/nthargout.m, scripts/general/num2str.m,
scripts/general/pol2cart.m, scripts/general/polyarea.m,
scripts/general/postpad.m, scripts/general/prepad.m,
scripts/general/private/__isequal__.m, scripts/general/private/__splinen__.m,
scripts/general/publish.m, scripts/general/quadgk.m, scripts/general/quadl.m,
scripts/general/quadv.m, scripts/general/rad2deg.m, scripts/general/randi.m,
scripts/general/rat.m, scripts/general/repmat.m, scripts/general/rot90.m,
scripts/general/rotdim.m, scripts/general/saveobj.m, scripts/general/shift.m,
scripts/general/shiftdim.m, scripts/general/sortrows.m,
scripts/general/sph2cart.m, scripts/general/structfun.m,
scripts/general/subsindex.m, scripts/general/trapz.m,
scripts/general/triplequad.m, scripts/general/validateattributes.m,
scripts/geometry/convhull.m, scripts/geometry/delaunay.m,
scripts/geometry/delaunayn.m, scripts/geometry/dsearch.m,
scripts/geometry/dsearchn.m, scripts/geometry/griddata.m,
scripts/geometry/griddata3.m, scripts/geometry/griddatan.m,
scripts/geometry/inpolygon.m, scripts/geometry/rectint.m,
scripts/geometry/tsearchn.m, scripts/geometry/voronoi.m,
scripts/geometry/voronoin.m, scripts/gui/dialog.m, scripts/gui/errordlg.m,
scripts/gui/guidata.m, scripts/gui/guihandles.m, scripts/gui/helpdlg.m,
scripts/gui/inputdlg.m, scripts/gui/listdlg.m, scripts/gui/msgbox.m,
scripts/gui/private/__file_filter__.m,
scripts/gui/private/__fltk_file_filter__.m,
scripts/gui/private/__get_funcname__.m, scripts/gui/private/__is_function__.m,
scripts/gui/private/__uigetdir_fltk__.m,
scripts/gui/private/__uigetfile_fltk__.m,
scripts/gui/private/__uiobject_split_args__.m,
scripts/gui/private/__uiputfile_fltk__.m, scripts/gui/questdlg.m,
scripts/gui/uibuttongroup.m, scripts/gui/uicontextmenu.m,
scripts/gui/uicontrol.m, scripts/gui/uigetdir.m, scripts/gui/uigetfile.m,
scripts/gui/uimenu.m, scripts/gui/uipanel.m, scripts/gui/uipushtool.m,
scripts/gui/uiputfile.m, scripts/gui/uiresume.m, scripts/gui/uitoggletool.m,
scripts/gui/uitoolbar.m, scripts/gui/uiwait.m, scripts/gui/waitbar.m,
scripts/gui/waitforbuttonpress.m, scripts/gui/warndlg.m,
scripts/help/__gripe_missing_component__.m, scripts/help/__makeinfo__.m,
scripts/help/__unimplemented__.m, scripts/help/ans.m, scripts/help/doc.m,
scripts/help/doc_cache_create.m, scripts/help/error_ids.m,
scripts/help/get_first_help_sentence.m, scripts/help/help.m,
scripts/help/lookfor.m, scripts/help/print_usage.m,
scripts/help/private/__additional_help_message__.m,
scripts/help/private/__strip_html_tags__.m, scripts/help/type.m,
scripts/help/warning_ids.m, scripts/help/which.m, scripts/image/autumn.m,
scripts/image/bone.m, scripts/image/brighten.m, scripts/image/cmpermute.m,
scripts/image/cmunique.m, scripts/image/colorcube.m, scripts/image/colormap.m,
scripts/image/contrast.m, scripts/image/cool.m, scripts/image/copper.m,
scripts/image/cubehelix.m, scripts/image/flag.m, scripts/image/gray.m,
scripts/image/gray2ind.m, scripts/image/hot.m, scripts/image/hsv.m,
scripts/image/hsv2rgb.m, scripts/image/im2double.m, scripts/image/image.m,
scripts/image/imagesc.m, scripts/image/imfinfo.m, scripts/image/imformats.m,
scripts/image/imread.m, scripts/image/imshow.m, scripts/image/imwrite.m,
scripts/image/ind2gray.m, scripts/image/ind2rgb.m, scripts/image/iscolormap.m,
scripts/image/jet.m, scripts/image/lines.m, scripts/image/ntsc2rgb.m,
scripts/image/ocean.m, scripts/image/pink.m, scripts/image/prism.m,
scripts/image/private/__imfinfo__.m, scripts/image/private/__imread__.m,
scripts/image/private/__imwrite__.m,
scripts/image/private/colorspace_conversion_input_check.m,
scripts/image/private/colorspace_conversion_revert.m,
scripts/image/private/imageIO.m, scripts/image/private/imwrite_filename.m,
scripts/image/private/ind2x.m, scripts/image/rainbow.m,
scripts/image/rgb2hsv.m, scripts/image/rgb2ind.m, scripts/image/rgb2ntsc.m,
scripts/image/rgbplot.m, scripts/image/spinmap.m, scripts/image/spring.m,
scripts/image/summer.m, scripts/image/viridis.m, scripts/image/white.m,
scripts/image/winter.m, scripts/io/beep.m, scripts/io/csvread.m,
scripts/io/csvwrite.m, scripts/io/dlmwrite.m, scripts/io/fileread.m,
scripts/io/importdata.m, scripts/io/is_valid_file_id.m, scripts/io/strread.m,
scripts/io/textread.m, scripts/java/javaArray.m, scripts/java/java_get.m,
scripts/java/java_set.m, scripts/java/javaaddpath.m, scripts/java/javachk.m,
scripts/java/javaclasspath.m, scripts/java/javamem.m,
scripts/java/javarmpath.m, scripts/java/org/octave/ClassHelper.java,
scripts/java/org/octave/Matrix.java,
scripts/java/org/octave/OctClassLoader.java,
scripts/java/org/octave/Octave.java,
scripts/java/org/octave/OctaveReference.java, scripts/java/usejava.m,
scripts/linear-algebra/bandwidth.m,
scripts/linear-algebra/commutation_matrix.m, scripts/linear-algebra/cond.m,
scripts/linear-algebra/condeig.m, scripts/linear-algebra/condest.m,
scripts/linear-algebra/cross.m, scripts/linear-algebra/duplication_matrix.m,
scripts/linear-algebra/expm.m, scripts/linear-algebra/housh.m,
scripts/linear-algebra/isbanded.m, scripts/linear-algebra/isdefinite.m,
scripts/linear-algebra/isdiag.m, scripts/linear-algebra/ishermitian.m,
scripts/linear-algebra/issymmetric.m, scripts/linear-algebra/istril.m,
scripts/linear-algebra/istriu.m, scripts/linear-algebra/krylov.m,
scripts/linear-algebra/linsolve.m, scripts/linear-algebra/logm.m,
scripts/linear-algebra/normest.m, scripts/linear-algebra/normest1.m,
scripts/linear-algebra/null.m, scripts/linear-algebra/orth.m,
scripts/linear-algebra/planerot.m, scripts/linear-algebra/qzhess.m,
scripts/linear-algebra/rank.m, scripts/linear-algebra/rref.m,
scripts/linear-algebra/subspace.m, scripts/linear-algebra/trace.m,
scripts/linear-algebra/vech.m, scripts/miscellaneous/bug_report.m,
scripts/miscellaneous/bunzip2.m, scripts/miscellaneous/cast.m,
scripts/miscellaneous/citation.m, scripts/miscellaneous/compare_versions.m,
scripts/miscellaneous/computer.m, scripts/miscellaneous/copyfile.m,
scripts/miscellaneous/debug.m, scripts/miscellaneous/delete.m,
scripts/miscellaneous/desktop.m, scripts/miscellaneous/dir.m,
scripts/miscellaneous/dos.m, scripts/miscellaneous/edit.m,
scripts/miscellaneous/fact.m, scripts/miscellaneous/fileattrib.m,
scripts/miscellaneous/fileparts.m, scripts/miscellaneous/fullfile.m,
scripts/miscellaneous/genvarname.m, scripts/miscellaneous/getappdata.m,
scripts/miscellaneous/getfield.m, scripts/miscellaneous/gunzip.m,
scripts/miscellaneous/info.m, scripts/miscellaneous/inputname.m,
scripts/miscellaneous/isappdata.m, scripts/miscellaneous/isdeployed.m,
scripts/miscellaneous/ismac.m, scripts/miscellaneous/ispc.m,
scripts/miscellaneous/isunix.m, scripts/miscellaneous/license.m,
scripts/miscellaneous/list_primes.m, scripts/miscellaneous/ls.m,
scripts/miscellaneous/ls_command.m, scripts/miscellaneous/menu.m,
scripts/miscellaneous/mex.m, scripts/miscellaneous/mexext.m,
scripts/miscellaneous/mkdir.m, scripts/miscellaneous/mkoctfile.m,
scripts/miscellaneous/movefile.m, scripts/miscellaneous/namelengthmax.m,
scripts/miscellaneous/news.m, scripts/miscellaneous/open.m,
scripts/miscellaneous/orderfields.m, scripts/miscellaneous/pack.m,
scripts/miscellaneous/parseparams.m, scripts/miscellaneous/perl.m,
scripts/miscellaneous/private/__w2mpth__.m,
scripts/miscellaneous/private/display_info_file.m,
scripts/miscellaneous/python.m, scripts/miscellaneous/recycle.m,
scripts/miscellaneous/rmappdata.m, scripts/miscellaneous/run.m,
scripts/miscellaneous/setappdata.m, scripts/miscellaneous/setfield.m,
scripts/miscellaneous/substruct.m, scripts/miscellaneous/swapbytes.m,
scripts/miscellaneous/symvar.m, scripts/miscellaneous/tar.m,
scripts/miscellaneous/tempdir.m, scripts/miscellaneous/tmpnam.m,
scripts/miscellaneous/unix.m, scripts/miscellaneous/unpack.m,
scripts/miscellaneous/untar.m, scripts/miscellaneous/unzip.m,
scripts/miscellaneous/ver.m, scripts/miscellaneous/version.m,
scripts/miscellaneous/what.m, scripts/miscellaneous/xor.m,
scripts/miscellaneous/zip.m, scripts/mk-pkg-add, scripts/mkdoc.pl,
scripts/ode/ode23.m, scripts/ode/ode45.m, scripts/ode/odeget.m,
scripts/ode/odeplot.m, scripts/ode/odeset.m, scripts/ode/private/AbsRel_norm.m,
scripts/ode/private/integrate_adaptive.m, scripts/ode/private/kahan.m,
scripts/ode/private/ode_event_handler.m, scripts/ode/private/odedefaults.m,
scripts/ode/private/odemergeopts.m, scripts/ode/private/runge_kutta_23.m,
scripts/ode/private/runge_kutta_45_dorpri.m,
scripts/ode/private/runge_kutta_interpolate.m,
scripts/ode/private/starting_stepsize.m, scripts/optimization/__all_opts__.m,
scripts/optimization/fminbnd.m, scripts/optimization/fminsearch.m,
scripts/optimization/fminunc.m, scripts/optimization/fsolve.m,
scripts/optimization/fzero.m, scripts/optimization/glpk.m,
scripts/optimization/lsqnonneg.m, scripts/optimization/optimget.m,
scripts/optimization/optimset.m, scripts/optimization/pqpnonneg.m,
scripts/optimization/private/__fdjac__.m, scripts/optimization/qp.m,
scripts/optimization/sqp.m, scripts/path/matlabroot.m, scripts/path/pathdef.m,
scripts/path/private/getsavepath.m, scripts/path/savepath.m, scripts/pkg/pkg.m,
scripts/pkg/private/build.m, scripts/pkg/private/configure_make.m,
scripts/pkg/private/default_prefix.m, scripts/pkg/private/describe.m,
scripts/pkg/private/dirempty.m, scripts/pkg/private/get_description.m,
scripts/pkg/private/get_forge_download.m, scripts/pkg/private/get_forge_pkg.m,
scripts/pkg/private/get_unsatisfied_deps.m, scripts/pkg/private/getarch.m,
scripts/pkg/private/getarchdir.m, scripts/pkg/private/install.m,
scripts/pkg/private/installed_packages.m,
scripts/pkg/private/list_forge_packages.m, scripts/pkg/private/load_packages.m,
scripts/pkg/private/load_packages_and_dependencies.m,
scripts/pkg/private/rebuild.m, scripts/pkg/private/save_order.m,
scripts/pkg/private/uninstall.m, scripts/pkg/private/unload_packages.m,
scripts/plot/appearance/__clabel__.m,
scripts/plot/appearance/__getlegenddata__.m,
scripts/plot/appearance/annotation.m, scripts/plot/appearance/axis.m,
scripts/plot/appearance/box.m, scripts/plot/appearance/caxis.m,
scripts/plot/appearance/clabel.m, scripts/plot/appearance/daspect.m,
scripts/plot/appearance/datetick.m, scripts/plot/appearance/diffuse.m,
scripts/plot/appearance/grid.m, scripts/plot/appearance/gtext.m,
scripts/plot/appearance/hidden.m, scripts/plot/appearance/legend.m,
scripts/plot/appearance/lighting.m, scripts/plot/appearance/material.m,
scripts/plot/appearance/orient.m, scripts/plot/appearance/pbaspect.m,
scripts/plot/appearance/private/__axis_label__.m,
scripts/plot/appearance/private/__axis_limits__.m,
scripts/plot/appearance/shading.m, scripts/plot/appearance/specular.m,
scripts/plot/appearance/text.m, scripts/plot/appearance/title.m,
scripts/plot/appearance/view.m, scripts/plot/appearance/whitebg.m,
scripts/plot/appearance/xlabel.m, scripts/plot/appearance/xlim.m,
scripts/plot/appearance/ylabel.m, scripts/plot/appearance/ylim.m,
scripts/plot/appearance/zlabel.m, scripts/plot/appearance/zlim.m,
scripts/plot/draw/area.m, scripts/plot/draw/bar.m, scripts/plot/draw/barh.m,
scripts/plot/draw/camlight.m, scripts/plot/draw/colorbar.m,
scripts/plot/draw/comet.m, scripts/plot/draw/comet3.m,
scripts/plot/draw/compass.m, scripts/plot/draw/contour.m,
scripts/plot/draw/contour3.m, scripts/plot/draw/contourc.m,
scripts/plot/draw/contourf.m, scripts/plot/draw/cylinder.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/ezmeshc.m,
scripts/plot/draw/ezplot.m, scripts/plot/draw/ezplot3.m,
scripts/plot/draw/ezpolar.m, scripts/plot/draw/ezsurf.m,
scripts/plot/draw/ezsurfc.m, scripts/plot/draw/feather.m,
scripts/plot/draw/fill.m, scripts/plot/draw/fplot.m, scripts/plot/draw/hist.m,
scripts/plot/draw/isocaps.m, scripts/plot/draw/isocolors.m,
scripts/plot/draw/isonormals.m, scripts/plot/draw/isosurface.m,
scripts/plot/draw/light.m, scripts/plot/draw/line.m,
scripts/plot/draw/loglog.m, scripts/plot/draw/loglogerr.m,
scripts/plot/draw/mesh.m, scripts/plot/draw/meshc.m, scripts/plot/draw/meshz.m,
scripts/plot/draw/pareto.m, scripts/plot/draw/patch.m,
scripts/plot/draw/pcolor.m, scripts/plot/draw/peaks.m, scripts/plot/draw/pie.m,
scripts/plot/draw/pie3.m, scripts/plot/draw/plot.m, scripts/plot/draw/plot3.m,
scripts/plot/draw/plotmatrix.m, scripts/plot/draw/plotyy.m,
scripts/plot/draw/polar.m, scripts/plot/draw/private/__add_datasource__.m,
scripts/plot/draw/private/__bar__.m,
scripts/plot/draw/private/__calc_isovalue_from_data__.m,
scripts/plot/draw/private/__contour__.m,
scripts/plot/draw/private/__errplot__.m,
scripts/plot/draw/private/__ezplot__.m,
scripts/plot/draw/private/__interp_cube__.m,
scripts/plot/draw/private/__line__.m,
scripts/plot/draw/private/__marching_cube__.m,
scripts/plot/draw/private/__patch__.m, scripts/plot/draw/private/__pie__.m,
scripts/plot/draw/private/__plt__.m, scripts/plot/draw/private/__quiver__.m,
scripts/plot/draw/private/__rotate_around_axis__.m,
scripts/plot/draw/private/__scatter__.m, scripts/plot/draw/private/__stem__.m,
scripts/plot/draw/private/__unite_shared_vertices__.m,
scripts/plot/draw/quiver.m, scripts/plot/draw/quiver3.m,
scripts/plot/draw/rectangle.m, scripts/plot/draw/reducepatch.m,
scripts/plot/draw/reducevolume.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/semilogxerr.m, scripts/plot/draw/semilogy.m,
scripts/plot/draw/semilogyerr.m, scripts/plot/draw/shrinkfaces.m,
scripts/plot/draw/slice.m, scripts/plot/draw/smooth3.m,
scripts/plot/draw/sombrero.m, scripts/plot/draw/sphere.m,
scripts/plot/draw/stairs.m, scripts/plot/draw/stem.m,
scripts/plot/draw/stem3.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/draw/waterfall.m,
scripts/plot/util/__actual_axis_position__.m,
scripts/plot/util/__default_plot_options__.m,
scripts/plot/util/__gnuplot_drawnow__.m,
scripts/plot/util/__next_line_color__.m,
scripts/plot/util/__next_line_style__.m, scripts/plot/util/__opengl_info__.m,
scripts/plot/util/__plt_get_axis_arg__.m, scripts/plot/util/__pltopt__.m,
scripts/plot/util/allchild.m, scripts/plot/util/ancestor.m,
scripts/plot/util/axes.m, scripts/plot/util/cla.m, scripts/plot/util/clf.m,
scripts/plot/util/close.m, scripts/plot/util/closereq.m,
scripts/plot/util/colstyle.m, scripts/plot/util/copyobj.m,
scripts/plot/util/figure.m, scripts/plot/util/findall.m,
scripts/plot/util/findfigs.m, scripts/plot/util/findobj.m,
scripts/plot/util/frame2im.m, scripts/plot/util/gca.m,
scripts/plot/util/gcbf.m, scripts/plot/util/gcbo.m, scripts/plot/util/gcf.m,
scripts/plot/util/gco.m, scripts/plot/util/ginput.m,
scripts/plot/util/gnuplot_binary.in.m, scripts/plot/util/graphics_toolkit.m,
scripts/plot/util/hdl2struct.m, scripts/plot/util/hggroup.m,
scripts/plot/util/hgload.m, scripts/plot/util/hgsave.m,
scripts/plot/util/hold.m, scripts/plot/util/im2frame.m,
scripts/plot/util/isaxes.m, scripts/plot/util/isfigure.m,
scripts/plot/util/ishghandle.m, scripts/plot/util/ishold.m,
scripts/plot/util/isprop.m, scripts/plot/util/linkaxes.m,
scripts/plot/util/linkprop.m, scripts/plot/util/meshgrid.m,
scripts/plot/util/ndgrid.m, scripts/plot/util/newplot.m,
scripts/plot/util/pan.m, scripts/plot/util/print.m, scripts/plot/util/printd.m,
scripts/plot/util/private/__add_default_menu__.m,
scripts/plot/util/private/__ghostscript__.m,
scripts/plot/util/private/__gnuplot_draw_axes__.m,
scripts/plot/util/private/__gnuplot_draw_figure__.m,
scripts/plot/util/private/__gnuplot_get_var__.m,
scripts/plot/util/private/__gnuplot_ginput__.m,
scripts/plot/util/private/__gnuplot_has_feature__.m,
scripts/plot/util/private/__gnuplot_has_terminal__.m,
scripts/plot/util/private/__gnuplot_open_stream__.m,
scripts/plot/util/private/__gnuplot_print__.m,
scripts/plot/util/private/__gnuplot_version__.m,
scripts/plot/util/private/__opengl_print__.m,
scripts/plot/util/private/__print_parse_opts__.m, scripts/plot/util/refresh.m,
scripts/plot/util/refreshdata.m, scripts/plot/util/rotate.m,
scripts/plot/util/rotate3d.m, scripts/plot/util/saveas.m,
scripts/plot/util/shg.m, scripts/plot/util/struct2hdl.m,
scripts/plot/util/subplot.m, scripts/plot/util/zoom.m,
scripts/polynomial/compan.m, scripts/polynomial/conv.m,
scripts/polynomial/deconv.m, scripts/polynomial/mkpp.m,
scripts/polynomial/mpoles.m, scripts/polynomial/padecoef.m,
scripts/polynomial/pchip.m, scripts/polynomial/poly.m,
scripts/polynomial/polyaffine.m, scripts/polynomial/polyder.m,
scripts/polynomial/polyeig.m, scripts/polynomial/polyfit.m,
scripts/polynomial/polygcd.m, scripts/polynomial/polyint.m,
scripts/polynomial/polyout.m, scripts/polynomial/polyreduce.m,
scripts/polynomial/polyval.m, scripts/polynomial/polyvalm.m,
scripts/polynomial/ppder.m, scripts/polynomial/ppint.m,
scripts/polynomial/ppjumps.m, scripts/polynomial/ppval.m,
scripts/polynomial/residue.m, scripts/polynomial/roots.m,
scripts/polynomial/spline.m, scripts/polynomial/splinefit.m,
scripts/polynomial/unmkpp.m, scripts/prefs/addpref.m, scripts/prefs/getpref.m,
scripts/prefs/ispref.m, scripts/prefs/prefdir.m, scripts/prefs/preferences.m,
scripts/prefs/private/loadprefs.m, scripts/prefs/private/prefsfile.m,
scripts/prefs/private/saveprefs.m, scripts/prefs/rmpref.m,
scripts/prefs/setpref.m, scripts/profiler/html/style.css,
scripts/profiler/profexplore.m, scripts/profiler/profexport.m,
scripts/profiler/profile.m, scripts/profiler/profshow.m,
scripts/set/intersect.m, scripts/set/ismember.m, scripts/set/powerset.m,
scripts/set/private/validsetargs.m, scripts/set/setdiff.m,
scripts/set/setxor.m, scripts/set/union.m, scripts/set/unique.m,
scripts/signal/arch_fit.m, scripts/signal/arch_rnd.m,
scripts/signal/arch_test.m, scripts/signal/arma_rnd.m,
scripts/signal/autoreg_matrix.m, scripts/signal/bartlett.m,
scripts/signal/blackman.m, scripts/signal/detrend.m, scripts/signal/diffpara.m,
scripts/signal/durbinlevinson.m, scripts/signal/fftconv.m,
scripts/signal/fftfilt.m, scripts/signal/fftshift.m, scripts/signal/filter2.m,
scripts/signal/fractdiff.m, scripts/signal/freqz.m,
scripts/signal/freqz_plot.m, scripts/signal/hamming.m,
scripts/signal/hanning.m, scripts/signal/hurst.m, scripts/signal/ifftshift.m,
scripts/signal/periodogram.m, scripts/signal/private/rectangle_lw.m,
scripts/signal/private/rectangle_sw.m, scripts/signal/private/triangle_lw.m,
scripts/signal/private/triangle_sw.m, scripts/signal/sinc.m,
scripts/signal/sinetone.m, scripts/signal/sinewave.m,
scripts/signal/spectral_adf.m, scripts/signal/spectral_xdf.m,
scripts/signal/spencer.m, scripts/signal/stft.m, scripts/signal/synthesis.m,
scripts/signal/unwrap.m, scripts/signal/yulewalker.m, scripts/sparse/bicg.m,
scripts/sparse/bicgstab.m, scripts/sparse/cgs.m, scripts/sparse/colperm.m,
scripts/sparse/eigs.m, scripts/sparse/etreeplot.m, scripts/sparse/gmres.m,
scripts/sparse/gplot.m, scripts/sparse/ichol.m, scripts/sparse/ilu.m,
scripts/sparse/nonzeros.m, scripts/sparse/pcg.m, scripts/sparse/pcr.m,
scripts/sparse/private/__sprand__.m, scripts/sparse/qmr.m,
scripts/sparse/spaugment.m, scripts/sparse/spconvert.m,
scripts/sparse/spdiags.m, scripts/sparse/speye.m, scripts/sparse/spfun.m,
scripts/sparse/spones.m, scripts/sparse/sprand.m, scripts/sparse/sprandn.m,
scripts/sparse/sprandsym.m, scripts/sparse/spstats.m, scripts/sparse/spy.m,
scripts/sparse/svds.m, scripts/sparse/treelayout.m, scripts/sparse/treeplot.m,
scripts/specfun/bessel.m, scripts/specfun/beta.m, scripts/specfun/betaln.m,
scripts/specfun/ellipke.m, scripts/specfun/expint.m, scripts/specfun/factor.m,
scripts/specfun/factorial.m, scripts/specfun/isprime.m, scripts/specfun/lcm.m,
scripts/specfun/legendre.m, scripts/specfun/nchoosek.m,
scripts/specfun/nthroot.m, scripts/specfun/perms.m, scripts/specfun/pow2.m,
scripts/specfun/primes.m, scripts/specfun/reallog.m, scripts/specfun/realpow.m,
scripts/specfun/realsqrt.m, scripts/special-matrix/gallery.m,
scripts/special-matrix/hadamard.m, scripts/special-matrix/hankel.m,
scripts/special-matrix/hilb.m, scripts/special-matrix/invhilb.m,
scripts/special-matrix/magic.m, scripts/special-matrix/pascal.m,
scripts/special-matrix/rosser.m, scripts/special-matrix/toeplitz.m,
scripts/special-matrix/vander.m, scripts/special-matrix/wilkinson.m,
scripts/startup/__finish__.m, scripts/statistics/base/center.m,
scripts/statistics/base/cloglog.m, scripts/statistics/base/corr.m,
scripts/statistics/base/cov.m, scripts/statistics/base/gls.m,
scripts/statistics/base/histc.m, scripts/statistics/base/iqr.m,
scripts/statistics/base/kendall.m, scripts/statistics/base/kurtosis.m,
scripts/statistics/base/logit.m, scripts/statistics/base/lscov.m,
scripts/statistics/base/mean.m, scripts/statistics/base/meansq.m,
scripts/statistics/base/median.m, scripts/statistics/base/mode.m,
scripts/statistics/base/moment.m, scripts/statistics/base/ols.m,
scripts/statistics/base/ppplot.m, scripts/statistics/base/prctile.m,
scripts/statistics/base/probit.m, scripts/statistics/base/qqplot.m,
scripts/statistics/base/quantile.m, scripts/statistics/base/range.m,
scripts/statistics/base/ranks.m, scripts/statistics/base/run_count.m,
scripts/statistics/base/runlength.m, scripts/statistics/base/skewness.m,
scripts/statistics/base/spearman.m, scripts/statistics/base/statistics.m,
scripts/statistics/base/std.m, scripts/statistics/base/table.m,
scripts/statistics/base/var.m, scripts/statistics/base/zscore.m,
scripts/statistics/distributions/betacdf.m,
scripts/statistics/distributions/betainv.m,
scripts/statistics/distributions/betapdf.m,
scripts/statistics/distributions/betarnd.m,
scripts/statistics/distributions/binocdf.m,
scripts/statistics/distributions/binoinv.m,
scripts/statistics/distributions/binopdf.m,
scripts/statistics/distributions/binornd.m,
scripts/statistics/distributions/cauchy_cdf.m,
scripts/statistics/distributions/cauchy_inv.m,
scripts/statistics/distributions/cauchy_pdf.m,
scripts/statistics/distributions/cauchy_rnd.m,
scripts/statistics/distributions/chi2cdf.m,
scripts/statistics/distributions/chi2inv.m,
scripts/statistics/distributions/chi2pdf.m,
scripts/statistics/distributions/chi2rnd.m,
scripts/statistics/distributions/discrete_cdf.m,
scripts/statistics/distributions/discrete_inv.m,
scripts/statistics/distributions/discrete_pdf.m,
scripts/statistics/distributions/discrete_rnd.m,
scripts/statistics/distributions/empirical_cdf.m,
scripts/statistics/distributions/empirical_inv.m,
scripts/statistics/distributions/empirical_pdf.m,
scripts/statistics/distributions/empirical_rnd.m,
scripts/statistics/distributions/expcdf.m,
scripts/statistics/distributions/expinv.m,
scripts/statistics/distributions/exppdf.m,
scripts/statistics/distributions/exprnd.m,
scripts/statistics/distributions/fcdf.m,
scripts/statistics/distributions/finv.m,
scripts/statistics/distributions/fpdf.m,
scripts/statistics/distributions/frnd.m,
scripts/statistics/distributions/gamcdf.m,
scripts/statistics/distributions/gaminv.m,
scripts/statistics/distributions/gampdf.m,
scripts/statistics/distributions/gamrnd.m,
scripts/statistics/distributions/geocdf.m,
scripts/statistics/distributions/geoinv.m,
scripts/statistics/distributions/geopdf.m,
scripts/statistics/distributions/geornd.m,
scripts/statistics/distributions/hygecdf.m,
scripts/statistics/distributions/hygeinv.m,
scripts/statistics/distributions/hygepdf.m,
scripts/statistics/distributions/hygernd.m,
scripts/statistics/distributions/kolmogorov_smirnov_cdf.m,
scripts/statistics/distributions/laplace_cdf.m,
scripts/statistics/distributions/laplace_inv.m,
scripts/statistics/distributions/laplace_pdf.m,
scripts/statistics/distributions/laplace_rnd.m,
scripts/statistics/distributions/logistic_cdf.m,
scripts/statistics/distributions/logistic_inv.m,
scripts/statistics/distributions/logistic_pdf.m,
scripts/statistics/distributions/logistic_rnd.m,
scripts/statistics/distributions/logncdf.m,
scripts/statistics/distributions/logninv.m,
scripts/statistics/distributions/lognpdf.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/normcdf.m,
scripts/statistics/distributions/norminv.m,
scripts/statistics/distributions/normpdf.m,
scripts/statistics/distributions/normrnd.m,
scripts/statistics/distributions/poisscdf.m,
scripts/statistics/distributions/poissinv.m,
scripts/statistics/distributions/poisspdf.m,
scripts/statistics/distributions/poissrnd.m,
scripts/statistics/distributions/stdnormal_cdf.m,
scripts/statistics/distributions/stdnormal_inv.m,
scripts/statistics/distributions/stdnormal_pdf.m,
scripts/statistics/distributions/stdnormal_rnd.m,
scripts/statistics/distributions/tcdf.m,
scripts/statistics/distributions/tinv.m,
scripts/statistics/distributions/tpdf.m,
scripts/statistics/distributions/trnd.m,
scripts/statistics/distributions/unidcdf.m,
scripts/statistics/distributions/unidinv.m,
scripts/statistics/distributions/unidpdf.m,
scripts/statistics/distributions/unidrnd.m,
scripts/statistics/distributions/unifcdf.m,
scripts/statistics/distributions/unifinv.m,
scripts/statistics/distributions/unifpdf.m,
scripts/statistics/distributions/unifrnd.m,
scripts/statistics/distributions/wblcdf.m,
scripts/statistics/distributions/wblinv.m,
scripts/statistics/distributions/wblpdf.m,
scripts/statistics/distributions/wblrnd.m,
scripts/statistics/distributions/wienrnd.m,
scripts/statistics/models/logistic_regression.m,
scripts/statistics/models/private/logistic_regression_derivatives.m,
scripts/statistics/models/private/logistic_regression_likelihood.m,
scripts/statistics/tests/anova.m, scripts/statistics/tests/bartlett_test.m,
scripts/statistics/tests/chisquare_test_homogeneity.m,
scripts/statistics/tests/chisquare_test_independence.m,
scripts/statistics/tests/cor_test.m,
scripts/statistics/tests/f_test_regression.m,
scripts/statistics/tests/hotelling_test.m,
scripts/statistics/tests/hotelling_test_2.m,
scripts/statistics/tests/kolmogorov_smirnov_test.m,
scripts/statistics/tests/kolmogorov_smirnov_test_2.m,
scripts/statistics/tests/kruskal_wallis_test.m,
scripts/statistics/tests/manova.m, scripts/statistics/tests/mcnemar_test.m,
scripts/statistics/tests/prop_test_2.m, scripts/statistics/tests/run_test.m,
scripts/statistics/tests/sign_test.m, scripts/statistics/tests/t_test.m,
scripts/statistics/tests/t_test_2.m,
scripts/statistics/tests/t_test_regression.m,
scripts/statistics/tests/u_test.m, scripts/statistics/tests/var_test.m,
scripts/statistics/tests/welch_test.m,
scripts/statistics/tests/wilcoxon_test.m, scripts/statistics/tests/z_test.m,
scripts/statistics/tests/z_test_2.m, scripts/strings/base2dec.m,
scripts/strings/bin2dec.m, scripts/strings/blanks.m, scripts/strings/cstrcat.m,
scripts/strings/deblank.m, scripts/strings/dec2base.m,
scripts/strings/dec2bin.m, scripts/strings/dec2hex.m,
scripts/strings/findstr.m, scripts/strings/hex2dec.m, scripts/strings/index.m,
scripts/strings/isletter.m, scripts/strings/isstrprop.m,
scripts/strings/mat2str.m, scripts/strings/ostrsplit.m,
scripts/strings/regexptranslate.m, scripts/strings/rindex.m,
scripts/strings/str2num.m, scripts/strings/strcat.m, scripts/strings/strchr.m,
scripts/strings/strjoin.m, scripts/strings/strjust.m,
scripts/strings/strmatch.m, scripts/strings/strsplit.m,
scripts/strings/strtok.m, scripts/strings/strtrim.m,
scripts/strings/strtrunc.m, scripts/strings/substr.m,
scripts/strings/untabify.m, scripts/strings/validatestring.m,
scripts/testfun/__have_feature__.m, scripts/testfun/__printf_assert__.m,
scripts/testfun/__prog_output_assert__.m, scripts/testfun/__run_test_suite__.m,
scripts/testfun/assert.m, scripts/testfun/demo.m, scripts/testfun/example.m,
scripts/testfun/fail.m, scripts/testfun/private/compare_plot_demos.m,
scripts/testfun/private/dump_demos.m,
scripts/testfun/private/html_compare_plot_demos.m, scripts/testfun/rundemos.m,
scripts/testfun/runtests.m, scripts/testfun/speed.m, scripts/testfun/test.m,
scripts/time/addtodate.m, scripts/time/asctime.m, scripts/time/calendar.m,
scripts/time/clock.m, scripts/time/ctime.m, scripts/time/date.m,
scripts/time/datenum.m, scripts/time/datestr.m, scripts/time/datevec.m,
scripts/time/eomday.m, scripts/time/etime.m, scripts/time/is_leap_year.m,
scripts/time/now.m, scripts/time/weekday.m, src/display-available.c,
src/display-available.h, src/main-cli.cc, src/main-gui.cc, src/main.in.cc,
src/mkoctfile.in.cc, src/octave-build-info.h, src/octave-build-info.in.cc,
src/octave-config.in.cc, src/shared-fcns.h, test/args.tst,
test/classdef/classdef.tst, test/classes/classes.tst, test/colormaps.tst,
test/complex.tst, test/ctor-vs-method/ctor-vs-method.tst, test/diag-perm.tst,
test/error.tst, test/eval-catch.tst,
test/fcn-handle-derived-resolution/fcn-handle-derived-resolution.tst,
test/fntests.m, test/for.tst, test/func.tst, test/global.tst, test/if.tst,
test/index.tst, test/io.tst, test/jit.tst, test/line-continue.tst,
test/logical-index.tst, test/nest/nest.tst, test/null-assign.tst,
test/parser.tst, test/prefer.tst, test/range.tst, test/recursion.tst,
test/return.tst, test/slice.tst, test/struct.tst, test/switch.tst,
test/system.tst, test/transpose.tst, test/try.tst, test/unwind.tst,
test/while.tst:
Use the same Copyright formatting throughout Octave.
author | Rik <rik@octave.org> |
---|---|
date | Sun, 13 Nov 2016 20:33:47 -0800 |
parents | 48c00363dc74 |
children | aa55d32100c9 |
line wrap: on
line source
/* Copyright (C) 2004-2016 David Bateman Copyright (C) 1998-2004 Andy Adler Copyright (C) 2010 VZLU Prague This file is part of Octave. Octave is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. Octave is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with Octave; see the file COPYING. If not, see <http://www.gnu.org/licenses/>. */ #if defined (HAVE_CONFIG_H) # include "config.h" #endif #include <cfloat> #include <iostream> #include <vector> #include "quit.h" #include "lo-ieee.h" #include "lo-mappers.h" #include "f77-fcn.h" #include "dRowVector.h" #include "lo-lapack-proto.h" #include "mx-m-cs.h" #include "mx-cs-m.h" #include "mx-cm-s.h" #include "mx-fcm-fs.h" #include "mx-s-cm.h" #include "mx-fs-fcm.h" #include "oct-locbuf.h" #include "dDiagMatrix.h" #include "CDiagMatrix.h" #include "CSparse.h" #include "boolSparse.h" #include "dSparse.h" #include "functor.h" #include "oct-spparms.h" #include "sparse-lu.h" #include "oct-sparse.h" #include "sparse-util.h" #include "sparse-chol.h" #include "sparse-qr.h" #include "Sparse-op-defs.h" #include "Sparse-diag-op-defs.h" #include "Sparse-perm-op-defs.h" // Define whether to use a basic QR solver or one that uses a Dulmange // Mendelsohn factorization to seperate the problem into under-determined, // well-determined and over-determined parts and solves them seperately #if ! defined (USE_QRSOLVE) # include "sparse-dmsolve.h" #endif SparseComplexMatrix::SparseComplexMatrix (const SparseMatrix& a) : MSparse<Complex> (a) { } SparseComplexMatrix::SparseComplexMatrix (const SparseBoolMatrix& a) : MSparse<Complex> (a.rows (), a.cols (), a.nnz ()) { octave_idx_type nc = cols (); octave_idx_type nz = a.nnz (); for (octave_idx_type i = 0; i < nc + 1; i++) cidx (i) = a.cidx (i); for (octave_idx_type i = 0; i < nz; i++) { data (i) = Complex (a.data (i)); ridx (i) = a.ridx (i); } } SparseComplexMatrix::SparseComplexMatrix (const ComplexDiagMatrix& a) : MSparse<Complex> (a.rows (), a.cols (), a.length ()) { octave_idx_type j = 0; octave_idx_type l = a.length (); for (octave_idx_type i = 0; i < l; i++) { cidx (i) = j; if (a(i, i) != 0.0) { data (j) = a(i, i); ridx (j) = i; j++; } } for (octave_idx_type i = l; i <= a.cols (); i++) cidx (i) = j; } bool SparseComplexMatrix::operator == (const SparseComplexMatrix& a) const { octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nz = nnz (); octave_idx_type nr_a = a.rows (); octave_idx_type nc_a = a.cols (); octave_idx_type nz_a = a.nnz (); if (nr != nr_a || nc != nc_a || nz != nz_a) return false; for (octave_idx_type i = 0; i < nc + 1; i++) if (cidx (i) != a.cidx (i)) return false; for (octave_idx_type i = 0; i < nz; i++) if (data (i) != a.data (i) || ridx (i) != a.ridx (i)) return false; return true; } bool SparseComplexMatrix::operator != (const SparseComplexMatrix& a) const { return !(*this == a); } bool SparseComplexMatrix::is_hermitian (void) const { octave_idx_type nr = rows (); octave_idx_type nc = cols (); if (nr == nc && nr > 0) { for (octave_idx_type j = 0; j < nc; j++) { for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { octave_idx_type ri = ridx (i); if (ri != j) { bool found = false; for (octave_idx_type k = cidx (ri); k < cidx (ri+1); k++) { if (ridx (k) == j) { if (data (i) == conj (data (k))) found = true; break; } } if (! found) return false; } } } return true; } return false; } static const Complex Complex_NaN_result (octave::numeric_limits<double>::NaN (), octave::numeric_limits<double>::NaN ()); SparseComplexMatrix SparseComplexMatrix::max (int dim) const { Array<octave_idx_type> dummy_idx; return max (dummy_idx, dim); } SparseComplexMatrix SparseComplexMatrix::max (Array<octave_idx_type>& idx_arg, int dim) const { SparseComplexMatrix result; dim_vector dv = dims (); octave_idx_type nr = dv(0); octave_idx_type nc = dv(1); if (dim >= dv.ndims ()) { idx_arg.resize (dim_vector (nr, nc), 0); return *this; } if (dim < 0) dim = dv.first_non_singleton (); if (dim == 0) { idx_arg.resize (dim_vector (nr == 0 ? 0 : 1, nc), 0); if (nr == 0 || nc == 0 || dim >= dv.ndims ()) return SparseComplexMatrix (nr == 0 ? 0 : 1, nc); octave_idx_type nel = 0; for (octave_idx_type j = 0; j < nc; j++) { Complex tmp_max; double abs_max = octave::numeric_limits<double>::NaN (); octave_idx_type idx_j = 0; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { if (ridx (i) != idx_j) break; else idx_j++; } if (idx_j != nr) { tmp_max = 0.; abs_max = 0.; } for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { Complex tmp = data (i); if (octave::math::isnan (tmp)) continue; double abs_tmp = std::abs (tmp); if (octave::math::isnan (abs_max) || abs_tmp > abs_max) { idx_j = ridx (i); tmp_max = tmp; abs_max = abs_tmp; } } idx_arg.elem (j) = octave::math::isnan (tmp_max) ? 0 : idx_j; if (abs_max != 0.) nel++; } result = SparseComplexMatrix (1, nc, nel); octave_idx_type ii = 0; result.xcidx (0) = 0; for (octave_idx_type j = 0; j < nc; j++) { Complex tmp = elem (idx_arg(j), j); if (tmp != 0.) { result.xdata (ii) = tmp; result.xridx (ii++) = 0; } result.xcidx (j+1) = ii; } } else { idx_arg.resize (dim_vector (nr, nc == 0 ? 0 : 1), 0); if (nr == 0 || nc == 0 || dim >= dv.ndims ()) return SparseComplexMatrix (nr, nc == 0 ? 0 : 1); OCTAVE_LOCAL_BUFFER (octave_idx_type, found, nr); for (octave_idx_type i = 0; i < nr; i++) found[i] = 0; for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) if (found[ridx (i)] == -j) found[ridx (i)] = -j - 1; for (octave_idx_type i = 0; i < nr; i++) if (found[i] > -nc && found[i] < 0) idx_arg.elem (i) = -found[i]; for (octave_idx_type j = 0; j < nc; j++) { for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { octave_idx_type ir = ridx (i); octave_idx_type ix = idx_arg.elem (ir); Complex tmp = data (i); if (octave::math::isnan (tmp)) continue; else if (ix == -1 || std::abs (tmp) > std::abs (elem (ir, ix))) idx_arg.elem (ir) = j; } } octave_idx_type nel = 0; for (octave_idx_type j = 0; j < nr; j++) if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.) nel++; result = SparseComplexMatrix (nr, 1, nel); octave_idx_type ii = 0; result.xcidx (0) = 0; result.xcidx (1) = nel; for (octave_idx_type j = 0; j < nr; j++) { if (idx_arg(j) == -1) { idx_arg(j) = 0; result.xdata (ii) = Complex_NaN_result; result.xridx (ii++) = j; } else { Complex tmp = elem (j, idx_arg(j)); if (tmp != 0.) { result.xdata (ii) = tmp; result.xridx (ii++) = j; } } } } return result; } SparseComplexMatrix SparseComplexMatrix::min (int dim) const { Array<octave_idx_type> dummy_idx; return min (dummy_idx, dim); } SparseComplexMatrix SparseComplexMatrix::min (Array<octave_idx_type>& idx_arg, int dim) const { SparseComplexMatrix result; dim_vector dv = dims (); octave_idx_type nr = dv(0); octave_idx_type nc = dv(1); if (dim >= dv.ndims ()) { idx_arg.resize (dim_vector (nr, nc), 0); return *this; } if (dim < 0) dim = dv.first_non_singleton (); if (dim == 0) { idx_arg.resize (dim_vector (nr == 0 ? 0 : 1, nc), 0); if (nr == 0 || nc == 0 || dim >= dv.ndims ()) return SparseComplexMatrix (nr == 0 ? 0 : 1, nc); octave_idx_type nel = 0; for (octave_idx_type j = 0; j < nc; j++) { Complex tmp_min; double abs_min = octave::numeric_limits<double>::NaN (); octave_idx_type idx_j = 0; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { if (ridx (i) != idx_j) break; else idx_j++; } if (idx_j != nr) { tmp_min = 0.; abs_min = 0.; } for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { Complex tmp = data (i); if (octave::math::isnan (tmp)) continue; double abs_tmp = std::abs (tmp); if (octave::math::isnan (abs_min) || abs_tmp < abs_min) { idx_j = ridx (i); tmp_min = tmp; abs_min = abs_tmp; } } idx_arg.elem (j) = octave::math::isnan (tmp_min) ? 0 : idx_j; if (abs_min != 0.) nel++; } result = SparseComplexMatrix (1, nc, nel); octave_idx_type ii = 0; result.xcidx (0) = 0; for (octave_idx_type j = 0; j < nc; j++) { Complex tmp = elem (idx_arg(j), j); if (tmp != 0.) { result.xdata (ii) = tmp; result.xridx (ii++) = 0; } result.xcidx (j+1) = ii; } } else { idx_arg.resize (dim_vector (nr, nc == 0 ? 0 : 1), 0); if (nr == 0 || nc == 0 || dim >= dv.ndims ()) return SparseComplexMatrix (nr, nc == 0 ? 0 : 1); OCTAVE_LOCAL_BUFFER (octave_idx_type, found, nr); for (octave_idx_type i = 0; i < nr; i++) found[i] = 0; for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) if (found[ridx (i)] == -j) found[ridx (i)] = -j - 1; for (octave_idx_type i = 0; i < nr; i++) if (found[i] > -nc && found[i] < 0) idx_arg.elem (i) = -found[i]; for (octave_idx_type j = 0; j < nc; j++) { for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { octave_idx_type ir = ridx (i); octave_idx_type ix = idx_arg.elem (ir); Complex tmp = data (i); if (octave::math::isnan (tmp)) continue; else if (ix == -1 || std::abs (tmp) < std::abs (elem (ir, ix))) idx_arg.elem (ir) = j; } } octave_idx_type nel = 0; for (octave_idx_type j = 0; j < nr; j++) if (idx_arg.elem (j) == -1 || elem (j, idx_arg.elem (j)) != 0.) nel++; result = SparseComplexMatrix (nr, 1, nel); octave_idx_type ii = 0; result.xcidx (0) = 0; result.xcidx (1) = nel; for (octave_idx_type j = 0; j < nr; j++) { if (idx_arg(j) == -1) { idx_arg(j) = 0; result.xdata (ii) = Complex_NaN_result; result.xridx (ii++) = j; } else { Complex tmp = elem (j, idx_arg(j)); if (tmp != 0.) { result.xdata (ii) = tmp; result.xridx (ii++) = j; } } } } return result; } /* %!assert (max (max (speye (65536) * 1i)), sparse (1i)) %!assert (min (min (speye (65536) * 1i)), sparse (0)) %!assert (size (max (sparse (8, 0), [], 1)), [1, 0]) %!assert (size (max (sparse (8, 0), [], 2)), [8, 0]) %!assert (size (max (sparse (0, 8), [], 1)), [0, 8]) %!assert (size (max (sparse (0, 8), [], 2)), [0, 1]) %!assert (size (min (sparse (8, 0), [], 1)), [1, 0]) %!assert (size (min (sparse (8, 0), [], 2)), [8, 0]) %!assert (size (min (sparse (0, 8), [], 1)), [0, 8]) %!assert (size (min (sparse (0, 8), [], 2)), [0, 1]) */ ComplexRowVector SparseComplexMatrix::row (octave_idx_type i) const { octave_idx_type nc = columns (); ComplexRowVector retval (nc, 0); for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type k = cidx (j); k < cidx (j+1); k++) { if (ridx (k) == i) { retval(j) = data (k); break; } } return retval; } ComplexColumnVector SparseComplexMatrix::column (octave_idx_type i) const { octave_idx_type nr = rows (); ComplexColumnVector retval (nr, 0); for (octave_idx_type k = cidx (i); k < cidx (i+1); k++) retval(ridx (k)) = data (k); return retval; } // destructive insert/delete/reorder operations SparseComplexMatrix& SparseComplexMatrix::insert (const SparseMatrix& a, octave_idx_type r, octave_idx_type c) { SparseComplexMatrix tmp (a); return insert (tmp /*a*/, r, c); } SparseComplexMatrix& SparseComplexMatrix::insert (const SparseComplexMatrix& a, octave_idx_type r, octave_idx_type c) { MSparse<Complex>::insert (a, r, c); return *this; } SparseComplexMatrix& SparseComplexMatrix::insert (const SparseMatrix& a, const Array<octave_idx_type>& indx) { SparseComplexMatrix tmp (a); return insert (tmp /*a*/, indx); } SparseComplexMatrix& SparseComplexMatrix::insert (const SparseComplexMatrix& a, const Array<octave_idx_type>& indx) { MSparse<Complex>::insert (a, indx); return *this; } SparseComplexMatrix SparseComplexMatrix::concat (const SparseComplexMatrix& rb, const Array<octave_idx_type>& ra_idx) { // Don't use numel to avoid all possiblity of an overflow if (rb.rows () > 0 && rb.cols () > 0) insert (rb, ra_idx(0), ra_idx(1)); return *this; } SparseComplexMatrix SparseComplexMatrix::concat (const SparseMatrix& rb, const Array<octave_idx_type>& ra_idx) { SparseComplexMatrix tmp (rb); if (rb.rows () > 0 && rb.cols () > 0) insert (tmp, ra_idx(0), ra_idx(1)); return *this; } ComplexMatrix SparseComplexMatrix::matrix_value (void) const { return Sparse<Complex>::array_value (); } SparseComplexMatrix SparseComplexMatrix::hermitian (void) const { octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nz = nnz (); SparseComplexMatrix retval (nc, nr, nz); for (octave_idx_type i = 0; i < nz; i++) retval.xcidx (ridx (i) + 1)++; // retval.xcidx[1:nr] holds the row degrees for rows 0:(nr-1) nz = 0; for (octave_idx_type i = 1; i <= nr; i++) { const octave_idx_type tmp = retval.xcidx (i); retval.xcidx (i) = nz; nz += tmp; } // retval.xcidx[1:nr] holds row entry *start* offsets for rows 0:(nr-1) for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type k = cidx (j); k < cidx (j+1); k++) { octave_idx_type q = retval.xcidx (ridx (k) + 1)++; retval.xridx (q) = j; retval.xdata (q) = conj (data (k)); } assert (nnz () == retval.xcidx (nr)); // retval.xcidx[1:nr] holds row entry *end* offsets for rows 0:(nr-1) // and retval.xcidx[0:(nr-1)] holds their row entry *start* offsets return retval; } SparseComplexMatrix conj (const SparseComplexMatrix& a) { octave_idx_type nr = a.rows (); octave_idx_type nc = a.cols (); octave_idx_type nz = a.nnz (); SparseComplexMatrix retval (nc, nr, nz); for (octave_idx_type i = 0; i < nc + 1; i++) retval.cidx (i) = a.cidx (i); for (octave_idx_type i = 0; i < nz; i++) { retval.data (i) = conj (a.data (i)); retval.ridx (i) = a.ridx (i); } return retval; } SparseComplexMatrix SparseComplexMatrix::inverse (void) const { octave_idx_type info; double rcond; MatrixType mattype (*this); return inverse (mattype, info, rcond, 0, 0); } SparseComplexMatrix SparseComplexMatrix::inverse (MatrixType& mattype) const { octave_idx_type info; double rcond; return inverse (mattype, info, rcond, 0, 0); } SparseComplexMatrix SparseComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info) const { double rcond; return inverse (mattype, info, rcond, 0, 0); } SparseComplexMatrix SparseComplexMatrix::dinverse (MatrixType &mattyp, octave_idx_type& info, double& rcond, const bool, const bool calccond) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); info = 0; if (nr == 0 || nc == 0 || nr != nc) (*current_liboctave_error_handler) ("inverse requires square matrix"); // Print spparms("spumoni") info if requested int typ = mattyp.type (); mattyp.info (); if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal) (*current_liboctave_error_handler) ("incorrect matrix type"); if (typ == MatrixType::Permuted_Diagonal) retval = transpose (); else retval = *this; // Force make_unique to be called Complex *v = retval.data (); if (calccond) { double dmax = 0.; double dmin = octave::numeric_limits<double>::Inf (); for (octave_idx_type i = 0; i < nr; i++) { double tmp = std::abs (v[i]); if (tmp > dmax) dmax = tmp; if (tmp < dmin) dmin = tmp; } rcond = dmin / dmax; } for (octave_idx_type i = 0; i < nr; i++) v[i] = 1.0 / v[i]; return retval; } SparseComplexMatrix SparseComplexMatrix::tinverse (MatrixType &mattyp, octave_idx_type& info, double& rcond, const bool, const bool calccond) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); info = 0; if (nr == 0 || nc == 0 || nr != nc) (*current_liboctave_error_handler) ("inverse requires square matrix"); // Print spparms("spumoni") info if requested int typ = mattyp.type (); mattyp.info (); if (typ != MatrixType::Upper && typ != MatrixType::Permuted_Upper && typ != MatrixType::Lower && typ != MatrixType::Permuted_Lower) (*current_liboctave_error_handler) ("incorrect matrix type"); double anorm = 0.; double ainvnorm = 0.; if (calccond) { // Calculate the 1-norm of matrix for rcond calculation for (octave_idx_type j = 0; j < nr; j++) { double atmp = 0.; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) atmp += std::abs (data (i)); if (atmp > anorm) anorm = atmp; } } if (typ == MatrixType::Upper || typ == MatrixType::Lower) { octave_idx_type nz = nnz (); octave_idx_type cx = 0; octave_idx_type nz2 = nz; retval = SparseComplexMatrix (nr, nc, nz2); for (octave_idx_type i = 0; i < nr; i++) { octave_quit (); // place the 1 in the identity position octave_idx_type cx_colstart = cx; if (cx == nz2) { nz2 *= 2; retval.change_capacity (nz2); } retval.xcidx (i) = cx; retval.xridx (cx) = i; retval.xdata (cx) = 1.0; cx++; // iterate accross columns of input matrix for (octave_idx_type j = i+1; j < nr; j++) { Complex v = 0.; // iterate to calculate sum octave_idx_type colXp = retval.xcidx (i); octave_idx_type colUp = cidx (j); octave_idx_type rpX, rpU; if (cidx (j) == cidx (j+1)) (*current_liboctave_error_handler) ("division by zero"); do { octave_quit (); rpX = retval.xridx (colXp); rpU = ridx (colUp); if (rpX < rpU) colXp++; else if (rpX > rpU) colUp++; else { v -= retval.xdata (colXp) * data (colUp); colXp++; colUp++; } } while (rpX < j && rpU < j && colXp < cx && colUp < nz); // get A(m,m) if (typ == MatrixType::Upper) colUp = cidx (j+1) - 1; else colUp = cidx (j); Complex pivot = data (colUp); if (pivot == 0. || ridx (colUp) != j) (*current_liboctave_error_handler) ("division by zero"); if (v != 0.) { if (cx == nz2) { nz2 *= 2; retval.change_capacity (nz2); } retval.xridx (cx) = j; retval.xdata (cx) = v / pivot; cx++; } } // get A(m,m) octave_idx_type colUp; if (typ == MatrixType::Upper) colUp = cidx (i+1) - 1; else colUp = cidx (i); Complex pivot = data (colUp); if (pivot == 0. || ridx (colUp) != i) (*current_liboctave_error_handler) ("division by zero"); if (pivot != 1.0) for (octave_idx_type j = cx_colstart; j < cx; j++) retval.xdata (j) /= pivot; } retval.xcidx (nr) = cx; retval.maybe_compress (); } else { octave_idx_type nz = nnz (); octave_idx_type cx = 0; octave_idx_type nz2 = nz; retval = SparseComplexMatrix (nr, nc, nz2); OCTAVE_LOCAL_BUFFER (Complex, work, nr); OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nr); octave_idx_type *perm = mattyp.triangular_perm (); if (typ == MatrixType::Permuted_Upper) { for (octave_idx_type i = 0; i < nr; i++) rperm[perm[i]] = i; } else { for (octave_idx_type i = 0; i < nr; i++) rperm[i] = perm[i]; for (octave_idx_type i = 0; i < nr; i++) perm[rperm[i]] = i; } for (octave_idx_type i = 0; i < nr; i++) { octave_quit (); octave_idx_type iidx = rperm[i]; for (octave_idx_type j = 0; j < nr; j++) work[j] = 0.; // place the 1 in the identity position work[iidx] = 1.0; // iterate accross columns of input matrix for (octave_idx_type j = iidx+1; j < nr; j++) { Complex v = 0.; octave_idx_type jidx = perm[j]; // iterate to calculate sum for (octave_idx_type k = cidx (jidx); k < cidx (jidx+1); k++) { octave_quit (); v -= work[ridx (k)] * data (k); } // get A(m,m) Complex pivot; if (typ == MatrixType::Permuted_Upper) pivot = data (cidx (jidx+1) - 1); else pivot = data (cidx (jidx)); if (pivot == 0.) (*current_liboctave_error_handler) ("division by zero"); work[j] = v / pivot; } // get A(m,m) octave_idx_type colUp; if (typ == MatrixType::Permuted_Upper) colUp = cidx (perm[iidx]+1) - 1; else colUp = cidx (perm[iidx]); Complex pivot = data (colUp); if (pivot == 0.) (*current_liboctave_error_handler) ("division by zero"); octave_idx_type new_cx = cx; for (octave_idx_type j = iidx; j < nr; j++) if (work[j] != 0.0) { new_cx++; if (pivot != 1.0) work[j] /= pivot; } if (cx < new_cx) { nz2 = (2*nz2 < new_cx ? new_cx : 2*nz2); retval.change_capacity (nz2); } retval.xcidx (i) = cx; for (octave_idx_type j = iidx; j < nr; j++) if (work[j] != 0.) { retval.xridx (cx) = j; retval.xdata (cx++) = work[j]; } } retval.xcidx (nr) = cx; retval.maybe_compress (); } if (calccond) { // Calculate the 1-norm of inverse matrix for rcond calculation for (octave_idx_type j = 0; j < nr; j++) { double atmp = 0.; for (octave_idx_type i = retval.cidx (j); i < retval.cidx (j+1); i++) atmp += std::abs (retval.data (i)); if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } return retval; } SparseComplexMatrix SparseComplexMatrix::inverse (MatrixType& mattype, octave_idx_type& info, double& rcond, bool, bool calc_cond) const { int typ = mattype.type (false); SparseComplexMatrix ret; if (typ == MatrixType::Unknown) typ = mattype.type (*this); if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal) ret = dinverse (mattype, info, rcond, true, calc_cond); else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) ret = tinverse (mattype, info, rcond, true, calc_cond).transpose (); else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) { MatrixType newtype = mattype.transpose (); ret = transpose ().tinverse (newtype, info, rcond, true, calc_cond); } else { if (mattype.is_hermitian ()) { MatrixType tmp_typ (MatrixType::Upper); octave::math::sparse_chol<SparseComplexMatrix> fact (*this, info, false); rcond = fact.rcond (); if (info == 0) { double rcond2; SparseMatrix Q = fact.Q (); SparseComplexMatrix InvL = fact.L ().transpose (). tinverse (tmp_typ, info, rcond2, true, false); ret = Q * InvL.hermitian () * InvL * Q.transpose (); } else { // Matrix is either singular or not positive definite mattype.mark_as_unsymmetric (); } } if (! mattype.is_hermitian ()) { octave_idx_type n = rows (); ColumnVector Qinit(n); for (octave_idx_type i = 0; i < n; i++) Qinit(i) = i; MatrixType tmp_typ (MatrixType::Upper); octave::math::sparse_lu<SparseComplexMatrix> fact (*this, Qinit, Matrix (), false, false); rcond = fact.rcond (); double rcond2; SparseComplexMatrix InvL = fact.L ().transpose (). tinverse (tmp_typ, info, rcond2, true, false); SparseComplexMatrix InvU = fact.U (). tinverse (tmp_typ, info, rcond2, true, false).transpose (); ret = fact.Pc ().transpose () * InvU * InvL * fact.Pr (); } } return ret; } ComplexDET SparseComplexMatrix::determinant (void) const { octave_idx_type info; double rcond; return determinant (info, rcond, 0); } ComplexDET SparseComplexMatrix::determinant (octave_idx_type& info) const { double rcond; return determinant (info, rcond, 0); } ComplexDET SparseComplexMatrix::determinant (octave_idx_type& err, double& rcond, bool) const { ComplexDET retval; #if defined (HAVE_UMFPACK) octave_idx_type nr = rows (); octave_idx_type nc = cols (); if (nr == 0 || nc == 0 || nr != nc) { retval = ComplexDET (1.0); } else { err = 0; // Setup the control parameters Matrix Control (UMFPACK_CONTROL, 1); double *control = Control.fortran_vec (); UMFPACK_ZNAME (defaults) (control); double tmp = octave_sparse_params::get_key ("spumoni"); if (! octave::math::isnan (tmp)) Control (UMFPACK_PRL) = tmp; tmp = octave_sparse_params::get_key ("piv_tol"); if (! octave::math::isnan (tmp)) { Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; Control (UMFPACK_PIVOT_TOLERANCE) = tmp; } // Set whether we are allowed to modify Q or not tmp = octave_sparse_params::get_key ("autoamd"); if (! octave::math::isnan (tmp)) Control (UMFPACK_FIXQ) = tmp; // Turn-off UMFPACK scaling for LU Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; UMFPACK_ZNAME (report_control) (control); const octave_idx_type *Ap = cidx (); const octave_idx_type *Ai = ridx (); const Complex *Ax = data (); UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai, reinterpret_cast<const double *> (Ax), 0, 1, control); void *Symbolic; Matrix Info (1, UMFPACK_INFO); double *info = Info.fortran_vec (); int status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai, reinterpret_cast<const double *> (Ax), 0, 0, &Symbolic, control, info); if (status < 0) { UMFPACK_ZNAME (report_status) (control, status); UMFPACK_ZNAME (report_info) (control, info); UMFPACK_ZNAME (free_symbolic) (&Symbolic); (*current_liboctave_error_handler) ("SparseComplexMatrix::determinant symbolic factorization failed"); } else { UMFPACK_ZNAME (report_symbolic) (Symbolic, control); void *Numeric; status = UMFPACK_ZNAME (numeric) (Ap, Ai, reinterpret_cast<const double *> (Ax), 0, Symbolic, &Numeric, control, info); UMFPACK_ZNAME (free_symbolic) (&Symbolic); rcond = Info (UMFPACK_RCOND); if (status < 0) { UMFPACK_ZNAME (report_status) (control, status); UMFPACK_ZNAME (report_info) (control, info); UMFPACK_ZNAME (free_numeric) (&Numeric); (*current_liboctave_error_handler) ("SparseComplexMatrix::determinant numeric factorization failed"); } else { UMFPACK_ZNAME (report_numeric) (Numeric, control); double c10[2], e10; status = UMFPACK_ZNAME (get_determinant) (c10, 0, &e10, Numeric, info); if (status < 0) { UMFPACK_ZNAME (report_status) (control, status); UMFPACK_ZNAME (report_info) (control, info); (*current_liboctave_error_handler) ("SparseComplexMatrix::determinant error calculating determinant"); } else retval = ComplexDET (Complex (c10[0], c10[1]), e10, 10); UMFPACK_ZNAME (free_numeric) (&Numeric); } } } #else octave_unused_parameter (err); octave_unused_parameter (rcond); (*current_liboctave_error_handler) ("support for UMFPACK was unavailable or disabled when liboctave was built"); #endif return retval; } ComplexMatrix SparseComplexMatrix::dsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler, bool calc_cond) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nm = (nc < nr ? nc : nr); err = 0; if (nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || nc == 0 || b.cols () == 0) retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); else { // Print spparms("spumoni") info if requested int typ = mattype.type (); mattype.info (); if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal) (*current_liboctave_error_handler) ("incorrect matrix type"); retval.resize (nc, b.cols (), Complex (0.,0.)); if (typ == MatrixType::Diagonal) for (octave_idx_type j = 0; j < b.cols (); j++) for (octave_idx_type i = 0; i < nm; i++) retval(i,j) = b(i,j) / data (i); else for (octave_idx_type j = 0; j < b.cols (); j++) for (octave_idx_type k = 0; k < nc; k++) for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) retval(k,j) = b(ridx (i),j) / data (i); if (calc_cond) { double dmax = 0.; double dmin = octave::numeric_limits<double>::Inf (); for (octave_idx_type i = 0; i < nm; i++) { double tmp = std::abs (data (i)); if (tmp > dmax) dmax = tmp; if (tmp < dmin) dmin = tmp; } rcond = dmin / dmax; } else rcond = 1.0; } return retval; } SparseComplexMatrix SparseComplexMatrix::dsolve (MatrixType &mattype, const SparseMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler, bool calc_cond) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nm = (nc < nr ? nc : nr); err = 0; if (nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || nc == 0 || b.cols () == 0) retval = SparseComplexMatrix (nc, b.cols ()); else { // Print spparms("spumoni") info if requested int typ = mattype.type (); mattype.info (); if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal) (*current_liboctave_error_handler) ("incorrect matrix type"); octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nnz (); retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx (0) = 0; octave_idx_type ii = 0; if (typ == MatrixType::Diagonal) for (octave_idx_type j = 0; j < b.cols (); j++) { for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) { if (b.ridx (i) >= nm) break; retval.xridx (ii) = b.ridx (i); retval.xdata (ii++) = b.data (i) / data (b.ridx (i)); } retval.xcidx (j+1) = ii; } else for (octave_idx_type j = 0; j < b.cols (); j++) { for (octave_idx_type l = 0; l < nc; l++) for (octave_idx_type i = cidx (l); i < cidx (l+1); i++) { bool found = false; octave_idx_type k; for (k = b.cidx (j); k < b.cidx (j+1); k++) if (ridx (i) == b.ridx (k)) { found = true; break; } if (found) { retval.xridx (ii) = l; retval.xdata (ii++) = b.data (k) / data (i); } } retval.xcidx (j+1) = ii; } if (calc_cond) { double dmax = 0.; double dmin = octave::numeric_limits<double>::Inf (); for (octave_idx_type i = 0; i < nm; i++) { double tmp = std::abs (data (i)); if (tmp > dmax) dmax = tmp; if (tmp < dmin) dmin = tmp; } rcond = dmin / dmax; } else rcond = 1.0; } return retval; } ComplexMatrix SparseComplexMatrix::dsolve (MatrixType &mattype, const ComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler, bool calc_cond) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nm = (nc < nr ? nc : nr); err = 0; if (nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || nc == 0 || b.cols () == 0) retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); else { // Print spparms("spumoni") info if requested int typ = mattype.type (); mattype.info (); if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal) (*current_liboctave_error_handler) ("incorrect matrix type"); retval.resize (nc, b.cols (), Complex (0.,0.)); if (typ == MatrixType::Diagonal) for (octave_idx_type j = 0; j < b.cols (); j++) for (octave_idx_type i = 0; i < nm; i++) retval(i,j) = b(i,j) / data (i); else for (octave_idx_type j = 0; j < b.cols (); j++) for (octave_idx_type k = 0; k < nc; k++) for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) retval(k,j) = b(ridx (i),j) / data (i); if (calc_cond) { double dmax = 0.; double dmin = octave::numeric_limits<double>::Inf (); for (octave_idx_type i = 0; i < nr; i++) { double tmp = std::abs (data (i)); if (tmp > dmax) dmax = tmp; if (tmp < dmin) dmin = tmp; } rcond = dmin / dmax; } else rcond = 1.0; } return retval; } SparseComplexMatrix SparseComplexMatrix::dsolve (MatrixType &mattype, const SparseComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler, bool calc_cond) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nm = (nc < nr ? nc : nr); err = 0; if (nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || nc == 0 || b.cols () == 0) retval = SparseComplexMatrix (nc, b.cols ()); else { // Print spparms("spumoni") info if requested int typ = mattype.type (); mattype.info (); if (typ != MatrixType::Diagonal && typ != MatrixType::Permuted_Diagonal) (*current_liboctave_error_handler) ("incorrect matrix type"); octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nnz (); retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx (0) = 0; octave_idx_type ii = 0; if (typ == MatrixType::Diagonal) for (octave_idx_type j = 0; j < b.cols (); j++) { for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) { if (b.ridx (i) >= nm) break; retval.xridx (ii) = b.ridx (i); retval.xdata (ii++) = b.data (i) / data (b.ridx (i)); } retval.xcidx (j+1) = ii; } else for (octave_idx_type j = 0; j < b.cols (); j++) { for (octave_idx_type l = 0; l < nc; l++) for (octave_idx_type i = cidx (l); i < cidx (l+1); i++) { bool found = false; octave_idx_type k; for (k = b.cidx (j); k < b.cidx (j+1); k++) if (ridx (i) == b.ridx (k)) { found = true; break; } if (found) { retval.xridx (ii) = l; retval.xdata (ii++) = b.data (k) / data (i); } } retval.xcidx (j+1) = ii; } if (calc_cond) { double dmax = 0.; double dmin = octave::numeric_limits<double>::Inf (); for (octave_idx_type i = 0; i < nm; i++) { double tmp = std::abs (data (i)); if (tmp > dmax) dmax = tmp; if (tmp < dmin) dmin = tmp; } rcond = dmin / dmax; } else rcond = 1.0; } return retval; } ComplexMatrix SparseComplexMatrix::utsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nm = (nc > nr ? nc : nr); err = 0; if (nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || nc == 0 || b.cols () == 0) retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); else { // Print spparms("spumoni") info if requested int typ = mattype.type (); mattype.info (); if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper) (*current_liboctave_error_handler) ("incorrect matrix type"); double anorm = 0.; double ainvnorm = 0.; octave_idx_type b_nc = b.cols (); rcond = 1.; if (calc_cond) { // Calculate the 1-norm of matrix for rcond calculation for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) atmp += std::abs (data (i)); if (atmp > anorm) anorm = atmp; } } if (typ == MatrixType::Permuted_Upper) { retval.resize (nc, b_nc); octave_idx_type *perm = mattype.triangular_perm (); OCTAVE_LOCAL_BUFFER (Complex, work, nm); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) work[i] = b(i,j); for (octave_idx_type i = nr; i < nc; i++) work[i] = 0.; for (octave_idx_type k = nc-1; k >= 0; k--) { octave_idx_type kidx = perm[k]; if (work[k] != 0.) { if (ridx (cidx (kidx+1)-1) != k || data (cidx (kidx+1)-1) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (cidx (kidx+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (kidx); i < cidx (kidx+1)-1; i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } for (octave_idx_type i = 0; i < nc; i++) retval(perm[i], j) = work[i]; } if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { octave_idx_type iidx = perm[k]; if (work[k] != 0.) { Complex tmp = work[k] / data (cidx (iidx+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (iidx); i < cidx (iidx+1)-1; i++) { octave_idx_type idx2 = ridx (i); work[idx2] = work[idx2] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = 0; i < j+1; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } else { OCTAVE_LOCAL_BUFFER (Complex, work, nm); retval.resize (nc, b_nc); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) work[i] = b(i,j); for (octave_idx_type i = nr; i < nc; i++) work[i] = 0.; for (octave_idx_type k = nc-1; k >= 0; k--) { if (work[k] != 0.) { if (ridx (cidx (k+1)-1) != k || data (cidx (k+1)-1) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (cidx (k+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } for (octave_idx_type i = 0; i < nc; i++) retval.xelem (i, j) = work[i]; } if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { if (work[k] != 0.) { Complex tmp = work[k] / data (cidx (k+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = 0; i < j+1; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } triangular_error: if (err != 0) { if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } return retval; } SparseComplexMatrix SparseComplexMatrix::utsolve (MatrixType &mattype, const SparseMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nm = (nc > nr ? nc : nr); err = 0; if (nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || nc == 0 || b.cols () == 0) retval = SparseComplexMatrix (nc, b.cols ()); else { // Print spparms("spumoni") info if requested int typ = mattype.type (); mattype.info (); if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper) (*current_liboctave_error_handler) ("incorrect matrix type"); double anorm = 0.; double ainvnorm = 0.; rcond = 1.; if (calc_cond) { // Calculate the 1-norm of matrix for rcond calculation for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) atmp += std::abs (data (i)); if (atmp > anorm) anorm = atmp; } } octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nnz (); retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx (0) = 0; octave_idx_type ii = 0; octave_idx_type x_nz = b_nz; if (typ == MatrixType::Permuted_Upper) { octave_idx_type *perm = mattype.triangular_perm (); OCTAVE_LOCAL_BUFFER (Complex, work, nm); OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc); for (octave_idx_type i = 0; i < nc; i++) rperm[perm[i]] = i; for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) work[b.ridx (i)] = b.data (i); for (octave_idx_type k = nc-1; k >= 0; k--) { octave_idx_type kidx = perm[k]; if (work[k] != 0.) { if (ridx (cidx (kidx+1)-1) != k || data (cidx (kidx+1)-1) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (cidx (kidx+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (kidx); i < cidx (kidx+1)-1; i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } // Count nonzeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) { // Resize the sparse matrix octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; retval.change_capacity (sz); x_nz = sz; } for (octave_idx_type i = 0; i < nc; i++) if (work[rperm[i]] != 0.) { retval.xridx (ii) = i; retval.xdata (ii++) = work[rperm[i]]; } retval.xcidx (j+1) = ii; } retval.maybe_compress (); if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { octave_idx_type iidx = perm[k]; if (work[k] != 0.) { Complex tmp = work[k] / data (cidx (iidx+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (iidx); i < cidx (iidx+1)-1; i++) { octave_idx_type idx2 = ridx (i); work[idx2] = work[idx2] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = 0; i < j+1; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } else { OCTAVE_LOCAL_BUFFER (Complex, work, nm); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) work[b.ridx (i)] = b.data (i); for (octave_idx_type k = nc-1; k >= 0; k--) { if (work[k] != 0.) { if (ridx (cidx (k+1)-1) != k || data (cidx (k+1)-1) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (cidx (k+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } // Count nonzeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) { // Resize the sparse matrix octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; retval.change_capacity (sz); x_nz = sz; } for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx (ii) = i; retval.xdata (ii++) = work[i]; } retval.xcidx (j+1) = ii; } retval.maybe_compress (); if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { if (work[k] != 0.) { Complex tmp = work[k] / data (cidx (k+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = 0; i < j+1; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } triangular_error: if (err != 0) { if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } return retval; } ComplexMatrix SparseComplexMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nm = (nc > nr ? nc : nr); err = 0; if (nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || nc == 0 || b.cols () == 0) retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); else { // Print spparms("spumoni") info if requested int typ = mattype.type (); mattype.info (); if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper) (*current_liboctave_error_handler) ("incorrect matrix type"); double anorm = 0.; double ainvnorm = 0.; octave_idx_type b_nc = b.cols (); rcond = 1.; if (calc_cond) { // Calculate the 1-norm of matrix for rcond calculation for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) atmp += std::abs (data (i)); if (atmp > anorm) anorm = atmp; } } if (typ == MatrixType::Permuted_Upper) { retval.resize (nc, b_nc); octave_idx_type *perm = mattype.triangular_perm (); OCTAVE_LOCAL_BUFFER (Complex, work, nm); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) work[i] = b(i,j); for (octave_idx_type i = nr; i < nc; i++) work[i] = 0.; for (octave_idx_type k = nc-1; k >= 0; k--) { octave_idx_type kidx = perm[k]; if (work[k] != 0.) { if (ridx (cidx (kidx+1)-1) != k || data (cidx (kidx+1)-1) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (cidx (kidx+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (kidx); i < cidx (kidx+1)-1; i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } for (octave_idx_type i = 0; i < nc; i++) retval(perm[i], j) = work[i]; } if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { octave_idx_type iidx = perm[k]; if (work[k] != 0.) { Complex tmp = work[k] / data (cidx (iidx+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (iidx); i < cidx (iidx+1)-1; i++) { octave_idx_type idx2 = ridx (i); work[idx2] = work[idx2] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = 0; i < j+1; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } else { OCTAVE_LOCAL_BUFFER (Complex, work, nm); retval.resize (nc, b_nc); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) work[i] = b(i,j); for (octave_idx_type i = nr; i < nc; i++) work[i] = 0.; for (octave_idx_type k = nc-1; k >= 0; k--) { if (work[k] != 0.) { if (ridx (cidx (k+1)-1) != k || data (cidx (k+1)-1) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (cidx (k+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } for (octave_idx_type i = 0; i < nc; i++) retval.xelem (i, j) = work[i]; } if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { if (work[k] != 0.) { Complex tmp = work[k] / data (cidx (k+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = 0; i < j+1; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } triangular_error: if (err != 0) { if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } return retval; } SparseComplexMatrix SparseComplexMatrix::utsolve (MatrixType &mattype, const SparseComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nm = (nc > nr ? nc : nr); err = 0; if (nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || nc == 0 || b.cols () == 0) retval = SparseComplexMatrix (nc, b.cols ()); else { // Print spparms("spumoni") info if requested int typ = mattype.type (); mattype.info (); if (typ != MatrixType::Permuted_Upper && typ != MatrixType::Upper) (*current_liboctave_error_handler) ("incorrect matrix type"); double anorm = 0.; double ainvnorm = 0.; rcond = 1.; if (calc_cond) { // Calculate the 1-norm of matrix for rcond calculation for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) atmp += std::abs (data (i)); if (atmp > anorm) anorm = atmp; } } octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nnz (); retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx (0) = 0; octave_idx_type ii = 0; octave_idx_type x_nz = b_nz; if (typ == MatrixType::Permuted_Upper) { octave_idx_type *perm = mattype.triangular_perm (); OCTAVE_LOCAL_BUFFER (Complex, work, nm); OCTAVE_LOCAL_BUFFER (octave_idx_type, rperm, nc); for (octave_idx_type i = 0; i < nc; i++) rperm[perm[i]] = i; for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) work[b.ridx (i)] = b.data (i); for (octave_idx_type k = nc-1; k >= 0; k--) { octave_idx_type kidx = perm[k]; if (work[k] != 0.) { if (ridx (cidx (kidx+1)-1) != k || data (cidx (kidx+1)-1) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (cidx (kidx+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (kidx); i < cidx (kidx+1)-1; i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } // Count nonzeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) { // Resize the sparse matrix octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; retval.change_capacity (sz); x_nz = sz; } for (octave_idx_type i = 0; i < nc; i++) if (work[rperm[i]] != 0.) { retval.xridx (ii) = i; retval.xdata (ii++) = work[rperm[i]]; } retval.xcidx (j+1) = ii; } retval.maybe_compress (); if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { octave_idx_type iidx = perm[k]; if (work[k] != 0.) { Complex tmp = work[k] / data (cidx (iidx+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (iidx); i < cidx (iidx+1)-1; i++) { octave_idx_type idx2 = ridx (i); work[idx2] = work[idx2] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = 0; i < j+1; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } else { OCTAVE_LOCAL_BUFFER (Complex, work, nm); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) work[b.ridx (i)] = b.data (i); for (octave_idx_type k = nr-1; k >= 0; k--) { if (work[k] != 0.) { if (ridx (cidx (k+1)-1) != k || data (cidx (k+1)-1) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (cidx (k+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } // Count nonzeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) { // Resize the sparse matrix octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; retval.change_capacity (sz); x_nz = sz; } for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx (ii) = i; retval.xdata (ii++) = work[i]; } retval.xcidx (j+1) = ii; } retval.maybe_compress (); if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = j; k >= 0; k--) { if (work[k] != 0.) { Complex tmp = work[k] / data (cidx (k+1)-1); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1)-1; i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = 0; i < j+1; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } triangular_error: if (err != 0) { if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } return retval; } ComplexMatrix SparseComplexMatrix::ltsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nm = (nc > nr ? nc : nr); err = 0; if (nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || nc == 0 || b.cols () == 0) retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); else { // Print spparms("spumoni") info if requested int typ = mattype.type (); mattype.info (); if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower) (*current_liboctave_error_handler) ("incorrect matrix type"); double anorm = 0.; double ainvnorm = 0.; octave_idx_type b_nc = b.cols (); rcond = 1.; if (calc_cond) { // Calculate the 1-norm of matrix for rcond calculation for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) atmp += std::abs (data (i)); if (atmp > anorm) anorm = atmp; } } if (typ == MatrixType::Permuted_Lower) { retval.resize (nc, b_nc); OCTAVE_LOCAL_BUFFER (Complex, work, nm); octave_idx_type *perm = mattype.triangular_perm (); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type i = 0; i < nr; i++) work[perm[i]] = b(i,j); for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { octave_idx_type minr = nr; octave_idx_type mini = 0; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) if (perm[ridx (i)] < minr) { minr = perm[ridx (i)]; mini = i; } if (minr != k || data (mini) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (mini); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) { if (i == mini) continue; octave_idx_type iidx = perm[ridx (i)]; work[iidx] = work[iidx] - tmp * data (i); } } } for (octave_idx_type i = 0; i < nc; i++) retval(i, j) = work[i]; } if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { octave_idx_type minr = nr; octave_idx_type mini = 0; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) if (perm[ridx (i)] < minr) { minr = perm[ridx (i)]; mini = i; } Complex tmp = work[k] / data (mini); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) { if (i == mini) continue; octave_idx_type iidx = perm[ridx (i)]; work[iidx] = work[iidx] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } else { OCTAVE_LOCAL_BUFFER (Complex, work, nm); retval.resize (nc, b_nc, 0.); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) work[i] = b(i,j); for (octave_idx_type i = nr; i < nc; i++) work[i] = 0.; for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { if (ridx (cidx (k)) != k || data (cidx (k)) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (cidx (k)); work[k] = tmp; for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } for (octave_idx_type i = 0; i < nc; i++) retval.xelem (i, j) = work[i]; } if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = j; k < nc; k++) { if (work[k] != 0.) { Complex tmp = work[k] / data (cidx (k)); work[k] = tmp; for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } triangular_error: if (err != 0) { if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } return retval; } SparseComplexMatrix SparseComplexMatrix::ltsolve (MatrixType &mattype, const SparseMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nm = (nc > nr ? nc : nr); err = 0; if (nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || nc == 0 || b.cols () == 0) retval = SparseComplexMatrix (nc, b.cols ()); else { // Print spparms("spumoni") info if requested int typ = mattype.type (); mattype.info (); if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower) (*current_liboctave_error_handler) ("incorrect matrix type"); double anorm = 0.; double ainvnorm = 0.; rcond = 1.; if (calc_cond) { // Calculate the 1-norm of matrix for rcond calculation for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) atmp += std::abs (data (i)); if (atmp > anorm) anorm = atmp; } } octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nnz (); retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx (0) = 0; octave_idx_type ii = 0; octave_idx_type x_nz = b_nz; if (typ == MatrixType::Permuted_Lower) { OCTAVE_LOCAL_BUFFER (Complex, work, nm); octave_idx_type *perm = mattype.triangular_perm (); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) work[perm[b.ridx (i)]] = b.data (i); for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { octave_idx_type minr = nr; octave_idx_type mini = 0; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) if (perm[ridx (i)] < minr) { minr = perm[ridx (i)]; mini = i; } if (minr != k || data (mini) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (mini); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) { if (i == mini) continue; octave_idx_type iidx = perm[ridx (i)]; work[iidx] = work[iidx] - tmp * data (i); } } } // Count nonzeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) { // Resize the sparse matrix octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; retval.change_capacity (sz); x_nz = sz; } for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx (ii) = i; retval.xdata (ii++) = work[i]; } retval.xcidx (j+1) = ii; } retval.maybe_compress (); if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { octave_idx_type minr = nr; octave_idx_type mini = 0; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) if (perm[ridx (i)] < minr) { minr = perm[ridx (i)]; mini = i; } Complex tmp = work[k] / data (mini); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) { if (i == mini) continue; octave_idx_type iidx = perm[ridx (i)]; work[iidx] = work[iidx] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } else { OCTAVE_LOCAL_BUFFER (Complex, work, nm); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) work[b.ridx (i)] = b.data (i); for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { if (ridx (cidx (k)) != k || data (cidx (k)) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (cidx (k)); work[k] = tmp; for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } // Count nonzeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) { // Resize the sparse matrix octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; retval.change_capacity (sz); x_nz = sz; } for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx (ii) = i; retval.xdata (ii++) = work[i]; } retval.xcidx (j+1) = ii; } retval.maybe_compress (); if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = j; k < nc; k++) { if (work[k] != 0.) { Complex tmp = work[k] / data (cidx (k)); work[k] = tmp; for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } triangular_error: if (err != 0) { if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } return retval; } ComplexMatrix SparseComplexMatrix::ltsolve (MatrixType &mattype, const ComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nm = (nc > nr ? nc : nr); err = 0; if (nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || nc == 0 || b.cols () == 0) retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); else { // Print spparms("spumoni") info if requested int typ = mattype.type (); mattype.info (); if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower) (*current_liboctave_error_handler) ("incorrect matrix type"); double anorm = 0.; double ainvnorm = 0.; octave_idx_type b_nc = b.cols (); rcond = 1.; if (calc_cond) { // Calculate the 1-norm of matrix for rcond calculation for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) atmp += std::abs (data (i)); if (atmp > anorm) anorm = atmp; } } if (typ == MatrixType::Permuted_Lower) { retval.resize (nc, b_nc); OCTAVE_LOCAL_BUFFER (Complex, work, nm); octave_idx_type *perm = mattype.triangular_perm (); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type i = 0; i < nr; i++) work[perm[i]] = b(i,j); for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { octave_idx_type minr = nr; octave_idx_type mini = 0; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) if (perm[ridx (i)] < minr) { minr = perm[ridx (i)]; mini = i; } if (minr != k || data (mini) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (mini); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) { if (i == mini) continue; octave_idx_type iidx = perm[ridx (i)]; work[iidx] = work[iidx] - tmp * data (i); } } } for (octave_idx_type i = 0; i < nc; i++) retval(i, j) = work[i]; } if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { octave_idx_type minr = nr; octave_idx_type mini = 0; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) if (perm[ridx (i)] < minr) { minr = perm[ridx (i)]; mini = i; } Complex tmp = work[k] / data (mini); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) { if (i == mini) continue; octave_idx_type iidx = perm[ridx (i)]; work[iidx] = work[iidx] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } else { OCTAVE_LOCAL_BUFFER (Complex, work, nm); retval.resize (nc, b_nc, 0.); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) work[i] = b(i,j); for (octave_idx_type i = nr; i < nc; i++) work[i] = 0.; for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { if (ridx (cidx (k)) != k || data (cidx (k)) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (cidx (k)); work[k] = tmp; for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } for (octave_idx_type i = 0; i < nc; i++) retval.xelem (i, j) = work[i]; } if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = j; k < nc; k++) { if (work[k] != 0.) { Complex tmp = work[k] / data (cidx (k)); work[k] = tmp; for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } triangular_error: if (err != 0) { if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } return retval; } SparseComplexMatrix SparseComplexMatrix::ltsolve (MatrixType &mattype, const SparseComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nm = (nc > nr ? nc : nr); err = 0; if (nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || nc == 0 || b.cols () == 0) retval = SparseComplexMatrix (nc, b.cols ()); else { // Print spparms("spumoni") info if requested int typ = mattype.type (); mattype.info (); if (typ != MatrixType::Permuted_Lower && typ != MatrixType::Lower) (*current_liboctave_error_handler) ("incorrect matrix type"); double anorm = 0.; double ainvnorm = 0.; rcond = 1.; if (calc_cond) { // Calculate the 1-norm of matrix for rcond calculation for (octave_idx_type j = 0; j < nc; j++) { double atmp = 0.; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) atmp += std::abs (data (i)); if (atmp > anorm) anorm = atmp; } } octave_idx_type b_nc = b.cols (); octave_idx_type b_nz = b.nnz (); retval = SparseComplexMatrix (nc, b_nc, b_nz); retval.xcidx (0) = 0; octave_idx_type ii = 0; octave_idx_type x_nz = b_nz; if (typ == MatrixType::Permuted_Lower) { OCTAVE_LOCAL_BUFFER (Complex, work, nm); octave_idx_type *perm = mattype.triangular_perm (); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) work[perm[b.ridx (i)]] = b.data (i); for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { octave_idx_type minr = nr; octave_idx_type mini = 0; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) if (perm[ridx (i)] < minr) { minr = perm[ridx (i)]; mini = i; } if (minr != k || data (mini) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (mini); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) { if (i == mini) continue; octave_idx_type iidx = perm[ridx (i)]; work[iidx] = work[iidx] - tmp * data (i); } } } // Count nonzeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) { // Resize the sparse matrix octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; retval.change_capacity (sz); x_nz = sz; } for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx (ii) = i; retval.xdata (ii++) = work[i]; } retval.xcidx (j+1) = ii; } retval.maybe_compress (); if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { octave_idx_type minr = nr; octave_idx_type mini = 0; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) if (perm[ridx (i)] < minr) { minr = perm[ridx (i)]; mini = i; } Complex tmp = work[k] / data (mini); work[k] = tmp; for (octave_idx_type i = cidx (k); i < cidx (k+1); i++) { if (i == mini) continue; octave_idx_type iidx = perm[ridx (i)]; work[iidx] = work[iidx] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } else { OCTAVE_LOCAL_BUFFER (Complex, work, nm); for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) work[b.ridx (i)] = b.data (i); for (octave_idx_type k = 0; k < nc; k++) { if (work[k] != 0.) { if (ridx (cidx (k)) != k || data (cidx (k)) == 0.) { err = -2; goto triangular_error; } Complex tmp = work[k] / data (cidx (k)); work[k] = tmp; for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } // Count nonzeros in work vector and adjust space in // retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) { // Resize the sparse matrix octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; retval.change_capacity (sz); x_nz = sz; } for (octave_idx_type i = 0; i < nc; i++) if (work[i] != 0.) { retval.xridx (ii) = i; retval.xdata (ii++) = work[i]; } retval.xcidx (j+1) = ii; } retval.maybe_compress (); if (calc_cond) { // Calculation of 1-norm of inv(*this) for (octave_idx_type i = 0; i < nm; i++) work[i] = 0.; for (octave_idx_type j = 0; j < nr; j++) { work[j] = 1.; for (octave_idx_type k = j; k < nc; k++) { if (work[k] != 0.) { Complex tmp = work[k] / data (cidx (k)); work[k] = tmp; for (octave_idx_type i = cidx (k)+1; i < cidx (k+1); i++) { octave_idx_type iidx = ridx (i); work[iidx] = work[iidx] - tmp * data (i); } } } double atmp = 0; for (octave_idx_type i = j; i < nc; i++) { atmp += std::abs (work[i]); work[i] = 0.; } if (atmp > ainvnorm) ainvnorm = atmp; } rcond = 1. / ainvnorm / anorm; } } triangular_error: if (err != 0) { if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } return retval; } ComplexMatrix SparseComplexMatrix::trisolve (MatrixType &mattype, const Matrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); err = 0; if (nr != nc || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || b.cols () == 0) retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); else if (calc_cond) (*current_liboctave_error_handler) ("calculation of condition number not implemented"); else { // Print spparms("spumoni") info if requested volatile int typ = mattype.type (); mattype.info (); if (typ == MatrixType::Tridiagonal_Hermitian) { OCTAVE_LOCAL_BUFFER (double, D, nr); OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); if (mattype.is_dense ()) { octave_idx_type ii = 0; for (octave_idx_type j = 0; j < nc-1; j++) { D[j] = octave::math::real (data (ii++)); DL[j] = data (ii); ii += 2; } D[nc-1] = octave::math::real (data (ii)); } else { D[0] = 0.; for (octave_idx_type i = 0; i < nr - 1; i++) { D[i+1] = 0.; DL[i] = 0.; } for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { if (ridx (i) == j) D[j] = octave::math::real (data (i)); else if (ridx (i) == j + 1) DL[j] = data (i); } } octave_idx_type b_nc = b.cols (); retval = ComplexMatrix (b); Complex *result = retval.fortran_vec (); F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, F77_DBLE_CMPLX_ARG (DL), F77_DBLE_CMPLX_ARG (result), b.rows (), err)); if (err != 0) { err = 0; mattype.mark_as_unsymmetric (); typ = MatrixType::Tridiagonal; } else rcond = 1.; } if (typ == MatrixType::Tridiagonal) { OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1); OCTAVE_LOCAL_BUFFER (Complex, D, nr); OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); if (mattype.is_dense ()) { octave_idx_type ii = 0; for (octave_idx_type j = 0; j < nc-1; j++) { D[j] = data (ii++); DL[j] = data (ii++); DU[j] = data (ii++); } D[nc-1] = data (ii); } else { D[0] = 0.; for (octave_idx_type i = 0; i < nr - 1; i++) { D[i+1] = 0.; DL[i] = 0.; DU[i] = 0.; } for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { if (ridx (i) == j) D[j] = data (i); else if (ridx (i) == j + 1) DL[j] = data (i); else if (ridx (i) == j - 1) DU[j-1] = data (i); } } octave_idx_type b_nc = b.cols (); retval = ComplexMatrix (b); Complex *result = retval.fortran_vec (); F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, F77_DBLE_CMPLX_ARG (DL), F77_DBLE_CMPLX_ARG (D), F77_DBLE_CMPLX_ARG (DU), F77_DBLE_CMPLX_ARG (result), b.rows (), err)); if (err != 0) { rcond = 0.; err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (); } else rcond = 1.; } else if (typ != MatrixType::Tridiagonal_Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; } SparseComplexMatrix SparseComplexMatrix::trisolve (MatrixType &mattype, const SparseMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); err = 0; if (nr != nc || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || b.cols () == 0) retval = SparseComplexMatrix (nc, b.cols ()); else if (calc_cond) (*current_liboctave_error_handler) ("calculation of condition number not implemented"); else { // Print spparms("spumoni") info if requested int typ = mattype.type (); mattype.info (); // Note can't treat symmetric case as there is no dpttrf function if (typ == MatrixType::Tridiagonal || typ == MatrixType::Tridiagonal_Hermitian) { OCTAVE_LOCAL_BUFFER (Complex, DU2, nr - 2); OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1); OCTAVE_LOCAL_BUFFER (Complex, D, nr); OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); Array<octave_idx_type> ipvt (dim_vector (nr, 1)); octave_idx_type *pipvt = ipvt.fortran_vec (); if (mattype.is_dense ()) { octave_idx_type ii = 0; for (octave_idx_type j = 0; j < nc-1; j++) { D[j] = data (ii++); DL[j] = data (ii++); DU[j] = data (ii++); } D[nc-1] = data (ii); } else { D[0] = 0.; for (octave_idx_type i = 0; i < nr - 1; i++) { D[i+1] = 0.; DL[i] = 0.; DU[i] = 0.; } for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { if (ridx (i) == j) D[j] = data (i); else if (ridx (i) == j + 1) DL[j] = data (i); else if (ridx (i) == j - 1) DU[j-1] = data (i); } } F77_XFCN (zgttrf, ZGTTRF, (nr, F77_DBLE_CMPLX_ARG (DL), F77_DBLE_CMPLX_ARG (D), F77_DBLE_CMPLX_ARG (DU), F77_DBLE_CMPLX_ARG (DU2), pipvt, err)); if (err != 0) { err = -2; rcond = 0.0; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (); } else { char job = 'N'; volatile octave_idx_type x_nz = b.nnz (); octave_idx_type b_nc = b.cols (); retval = SparseComplexMatrix (nr, b_nc, x_nz); retval.xcidx (0) = 0; volatile octave_idx_type ii = 0; rcond = 1.0; OCTAVE_LOCAL_BUFFER (Complex, work, nr); for (volatile octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) work[i] = 0.; for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) work[b.ridx (i)] = b.data (i); F77_XFCN (zgttrs, ZGTTRS, (F77_CONST_CHAR_ARG2 (&job, 1), nr, 1, F77_DBLE_CMPLX_ARG (DL), F77_DBLE_CMPLX_ARG (D), F77_DBLE_CMPLX_ARG (DU), F77_DBLE_CMPLX_ARG (DU2), pipvt, F77_DBLE_CMPLX_ARG (work), b.rows (), err F77_CHAR_ARG_LEN (1))); // Count nonzeros in work vector and adjust // space in retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nr; i++) if (work[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) { // Resize the sparse matrix octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; retval.change_capacity (sz); x_nz = sz; } for (octave_idx_type i = 0; i < nr; i++) if (work[i] != 0.) { retval.xridx (ii) = i; retval.xdata (ii++) = work[i]; } retval.xcidx (j+1) = ii; } retval.maybe_compress (); } } else if (typ != MatrixType::Tridiagonal_Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; } ComplexMatrix SparseComplexMatrix::trisolve (MatrixType &mattype, const ComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); err = 0; if (nr != nc || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || b.cols () == 0) retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); else if (calc_cond) (*current_liboctave_error_handler) ("calculation of condition number not implemented"); else { // Print spparms("spumoni") info if requested volatile int typ = mattype.type (); mattype.info (); if (typ == MatrixType::Tridiagonal_Hermitian) { OCTAVE_LOCAL_BUFFER (double, D, nr); OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); if (mattype.is_dense ()) { octave_idx_type ii = 0; for (octave_idx_type j = 0; j < nc-1; j++) { D[j] = octave::math::real (data (ii++)); DL[j] = data (ii); ii += 2; } D[nc-1] = octave::math::real (data (ii)); } else { D[0] = 0.; for (octave_idx_type i = 0; i < nr - 1; i++) { D[i+1] = 0.; DL[i] = 0.; } for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { if (ridx (i) == j) D[j] = octave::math::real (data (i)); else if (ridx (i) == j + 1) DL[j] = data (i); } } octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); rcond = 1.; retval = ComplexMatrix (b); Complex *result = retval.fortran_vec (); F77_XFCN (zptsv, ZPTSV, (nr, b_nc, D, F77_DBLE_CMPLX_ARG (DL), F77_DBLE_CMPLX_ARG (result), b_nr, err)); if (err != 0) { err = 0; mattype.mark_as_unsymmetric (); typ = MatrixType::Tridiagonal; } } if (typ == MatrixType::Tridiagonal) { OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1); OCTAVE_LOCAL_BUFFER (Complex, D, nr); OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); if (mattype.is_dense ()) { octave_idx_type ii = 0; for (octave_idx_type j = 0; j < nc-1; j++) { D[j] = data (ii++); DL[j] = data (ii++); DU[j] = data (ii++); } D[nc-1] = data (ii); } else { D[0] = 0.; for (octave_idx_type i = 0; i < nr - 1; i++) { D[i+1] = 0.; DL[i] = 0.; DU[i] = 0.; } for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { if (ridx (i) == j) D[j] = data (i); else if (ridx (i) == j + 1) DL[j] = data (i); else if (ridx (i) == j - 1) DU[j-1] = data (i); } } octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); rcond = 1.; retval = ComplexMatrix (b); Complex *result = retval.fortran_vec (); F77_XFCN (zgtsv, ZGTSV, (nr, b_nc, F77_DBLE_CMPLX_ARG (DL), F77_DBLE_CMPLX_ARG (D), F77_DBLE_CMPLX_ARG (DU), F77_DBLE_CMPLX_ARG (result), b_nr, err)); if (err != 0) { rcond = 0.; err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (); } } else if (typ != MatrixType::Tridiagonal_Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; } SparseComplexMatrix SparseComplexMatrix::trisolve (MatrixType &mattype, const SparseComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); err = 0; if (nr != nc || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || b.cols () == 0) retval = SparseComplexMatrix (nc, b.cols ()); else if (calc_cond) (*current_liboctave_error_handler) ("calculation of condition number not implemented"); else { // Print spparms("spumoni") info if requested int typ = mattype.type (); mattype.info (); // Note can't treat symmetric case as there is no dpttrf function if (typ == MatrixType::Tridiagonal || typ == MatrixType::Tridiagonal_Hermitian) { OCTAVE_LOCAL_BUFFER (Complex, DU2, nr - 2); OCTAVE_LOCAL_BUFFER (Complex, DU, nr - 1); OCTAVE_LOCAL_BUFFER (Complex, D, nr); OCTAVE_LOCAL_BUFFER (Complex, DL, nr - 1); Array<octave_idx_type> ipvt (dim_vector (nr, 1)); octave_idx_type *pipvt = ipvt.fortran_vec (); if (mattype.is_dense ()) { octave_idx_type ii = 0; for (octave_idx_type j = 0; j < nc-1; j++) { D[j] = data (ii++); DL[j] = data (ii++); DU[j] = data (ii++); } D[nc-1] = data (ii); } else { D[0] = 0.; for (octave_idx_type i = 0; i < nr - 1; i++) { D[i+1] = 0.; DL[i] = 0.; DU[i] = 0.; } for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { if (ridx (i) == j) D[j] = data (i); else if (ridx (i) == j + 1) DL[j] = data (i); else if (ridx (i) == j - 1) DU[j-1] = data (i); } } F77_XFCN (zgttrf, ZGTTRF, (nr, F77_DBLE_CMPLX_ARG (DL), F77_DBLE_CMPLX_ARG (D), F77_DBLE_CMPLX_ARG (DU), F77_DBLE_CMPLX_ARG (DU2), pipvt, err)); if (err != 0) { rcond = 0.0; err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (); } else { rcond = 1.; char job = 'N'; octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); // Take a first guess that the number of nonzero terms // will be as many as in b volatile octave_idx_type x_nz = b.nnz (); volatile octave_idx_type ii = 0; retval = SparseComplexMatrix (b_nr, b_nc, x_nz); retval.xcidx (0) = 0; for (volatile octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < b_nr; i++) Bx[i] = b(i,j); F77_XFCN (zgttrs, ZGTTRS, (F77_CONST_CHAR_ARG2 (&job, 1), nr, 1, F77_DBLE_CMPLX_ARG (DL), F77_DBLE_CMPLX_ARG (D), F77_DBLE_CMPLX_ARG (DU), F77_DBLE_CMPLX_ARG (DU2), pipvt, F77_DBLE_CMPLX_ARG (Bx), b_nr, err F77_CHAR_ARG_LEN (1))); if (err != 0) { // FIXME: This should probably be a warning so that // error value can be passed back. (*current_liboctave_error_handler) ("SparseComplexMatrix::solve solve failed"); err = -1; break; } // Count nonzeros in work vector and adjust // space in retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nr; i++) if (Bx[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) { // Resize the sparse matrix octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; retval.change_capacity (sz); x_nz = sz; } for (octave_idx_type i = 0; i < nr; i++) if (Bx[i] != 0.) { retval.xridx (ii) = i; retval.xdata (ii++) = Bx[i]; } retval.xcidx (j+1) = ii; } retval.maybe_compress (); } } else if (typ != MatrixType::Tridiagonal_Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; } ComplexMatrix SparseComplexMatrix::bsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); err = 0; if (nr != nc || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || b.cols () == 0) retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); else { // Print spparms("spumoni") info if requested volatile int typ = mattype.type (); mattype.info (); if (typ == MatrixType::Banded_Hermitian) { octave_idx_type n_lower = mattype.nlower (); octave_idx_type ldm = n_lower + 1; ComplexMatrix m_band (ldm, nc); Complex *tmp_data = m_band.fortran_vec (); if (! mattype.is_dense ()) { octave_idx_type ii = 0; for (octave_idx_type j = 0; j < ldm; j++) for (octave_idx_type i = 0; i < nc; i++) tmp_data[ii++] = 0.; } for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { octave_idx_type ri = ridx (i); if (ri >= j) m_band(ri - j, j) = data (i); } // Calculate the norm of the matrix, for later use. double anorm; if (calc_cond) anorm = m_band.abs ().sum ().row (0).max (); char job = 'L'; F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, F77_DBLE_CMPLX_ARG (tmp_data), ldm, err F77_CHAR_ARG_LEN (1))); if (err != 0) { rcond = 0.0; // Matrix is not positive definite!! Fall through to // unsymmetric banded solver. mattype.mark_as_unsymmetric (); typ = MatrixType::Banded; err = 0; } else { if (calc_cond) { Array<Complex> z (dim_vector (2 * nr, 1)); Complex *pz = z.fortran_vec (); Array<double> iz (dim_vector (nr, 1)); double *piz = iz.fortran_vec (); F77_XFCN (zpbcon, ZPBCON, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, F77_DBLE_CMPLX_ARG (tmp_data), ldm, anorm, rcond, F77_DBLE_CMPLX_ARG (pz), piz, err F77_CHAR_ARG_LEN (1))); if (err != 0) err = -2; volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } else rcond = 1.0; if (err == 0) { retval = ComplexMatrix (b); Complex *result = retval.fortran_vec (); octave_idx_type b_nc = b.cols (); F77_XFCN (zpbtrs, ZPBTRS, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, b_nc, F77_DBLE_CMPLX_ARG (tmp_data), ldm, F77_DBLE_CMPLX_ARG (result), b.rows (), err F77_CHAR_ARG_LEN (1))); if (err != 0) { // FIXME: Probably should be a warning. (*current_liboctave_error_handler) ("SparseMatrix::solve solve failed"); err = -1; } } } } if (typ == MatrixType::Banded) { // Create the storage for the banded form of the sparse matrix octave_idx_type n_upper = mattype.nupper (); octave_idx_type n_lower = mattype.nlower (); octave_idx_type ldm = n_upper + 2 * n_lower + 1; ComplexMatrix m_band (ldm, nc); Complex *tmp_data = m_band.fortran_vec (); if (! mattype.is_dense ()) { octave_idx_type ii = 0; for (octave_idx_type j = 0; j < ldm; j++) for (octave_idx_type i = 0; i < nc; i++) tmp_data[ii++] = 0.; } for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) m_band(ridx (i) - j + n_lower + n_upper, j) = data (i); // Calculate the norm of the matrix, for later use. double anorm = 0.0; if (calc_cond) { for (octave_idx_type j = 0; j < nr; j++) { double atmp = 0.; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) atmp += std::abs (data (i)); if (atmp > anorm) anorm = atmp; } } Array<octave_idx_type> ipvt (dim_vector (nr, 1)); octave_idx_type *pipvt = ipvt.fortran_vec (); F77_XFCN (zgbtrf, ZGBTRF, (nr, nc, n_lower, n_upper, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt, err)); // Throw-away extra info LAPACK gives so as to not // change output. if (err != 0) { rcond = 0.0; err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (); } else { if (calc_cond) { char job = '1'; Array<Complex> z (dim_vector (2 * nr, 1)); Complex *pz = z.fortran_vec (); Array<double> iz (dim_vector (nr, 1)); double *piz = iz.fortran_vec (); F77_XFCN (zgbcon, ZGBCON, (F77_CONST_CHAR_ARG2 (&job, 1), nc, n_lower, n_upper, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt, anorm, rcond, F77_DBLE_CMPLX_ARG (pz), piz, err F77_CHAR_ARG_LEN (1))); if (err != 0) err = -2; volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } else rcond = 1.; if (err == 0) { retval = ComplexMatrix (b); Complex *result = retval.fortran_vec (); octave_idx_type b_nc = b.cols (); char job = 'N'; F77_XFCN (zgbtrs, ZGBTRS, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, n_upper, b_nc, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt, F77_DBLE_CMPLX_ARG (result), b.rows (), err F77_CHAR_ARG_LEN (1))); } } } else if (typ != MatrixType::Banded_Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; } SparseComplexMatrix SparseComplexMatrix::bsolve (MatrixType &mattype, const SparseMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); err = 0; if (nr != nc || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || b.cols () == 0) retval = SparseComplexMatrix (nc, b.cols ()); else { // Print spparms("spumoni") info if requested volatile int typ = mattype.type (); mattype.info (); if (typ == MatrixType::Banded_Hermitian) { octave_idx_type n_lower = mattype.nlower (); octave_idx_type ldm = n_lower + 1; ComplexMatrix m_band (ldm, nc); Complex *tmp_data = m_band.fortran_vec (); if (! mattype.is_dense ()) { octave_idx_type ii = 0; for (octave_idx_type j = 0; j < ldm; j++) for (octave_idx_type i = 0; i < nc; i++) tmp_data[ii++] = 0.; } for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { octave_idx_type ri = ridx (i); if (ri >= j) m_band(ri - j, j) = data (i); } // Calculate the norm of the matrix, for later use. double anorm; if (calc_cond) anorm = m_band.abs ().sum ().row (0).max (); char job = 'L'; F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, F77_DBLE_CMPLX_ARG (tmp_data), ldm, err F77_CHAR_ARG_LEN (1))); if (err != 0) { rcond = 0.0; mattype.mark_as_unsymmetric (); typ = MatrixType::Banded; err = 0; } else { if (calc_cond) { Array<Complex> z (dim_vector (2 * nr, 1)); Complex *pz = z.fortran_vec (); Array<double> iz (dim_vector (nr, 1)); double *piz = iz.fortran_vec (); F77_XFCN (zpbcon, ZPBCON, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, F77_DBLE_CMPLX_ARG (tmp_data), ldm, anorm, rcond, F77_DBLE_CMPLX_ARG (pz), piz, err F77_CHAR_ARG_LEN (1))); if (err != 0) err = -2; volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } else rcond = 1.0; if (err == 0) { octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); // Take a first guess that the number of nonzero terms // will be as many as in b volatile octave_idx_type x_nz = b.nnz (); volatile octave_idx_type ii = 0; retval = SparseComplexMatrix (b_nr, b_nc, x_nz); retval.xcidx (0) = 0; for (volatile octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < b_nr; i++) Bx[i] = b.elem (i, j); F77_XFCN (zpbtrs, ZPBTRS, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, 1, F77_DBLE_CMPLX_ARG (tmp_data), ldm, F77_DBLE_CMPLX_ARG (Bx), b_nr, err F77_CHAR_ARG_LEN (1))); if (err != 0) { // FIXME: Probably should be a warning. (*current_liboctave_error_handler) ("SparseComplexMatrix::solve solve failed"); err = -1; break; } for (octave_idx_type i = 0; i < b_nr; i++) { Complex tmp = Bx[i]; if (tmp != 0.0) { if (ii == x_nz) { // Resize the sparse matrix octave_idx_type sz = x_nz * (b_nc - j) / b_nc; sz = (sz > 10 ? sz : 10) + x_nz; retval.change_capacity (sz); x_nz = sz; } retval.xdata (ii) = tmp; retval.xridx (ii++) = i; } } retval.xcidx (j+1) = ii; } retval.maybe_compress (); } } } if (typ == MatrixType::Banded) { // Create the storage for the banded form of the sparse matrix octave_idx_type n_upper = mattype.nupper (); octave_idx_type n_lower = mattype.nlower (); octave_idx_type ldm = n_upper + 2 * n_lower + 1; ComplexMatrix m_band (ldm, nc); Complex *tmp_data = m_band.fortran_vec (); if (! mattype.is_dense ()) { octave_idx_type ii = 0; for (octave_idx_type j = 0; j < ldm; j++) for (octave_idx_type i = 0; i < nc; i++) tmp_data[ii++] = 0.; } for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) m_band(ridx (i) - j + n_lower + n_upper, j) = data (i); // Calculate the norm of the matrix, for later use. double anorm = 0.0; if (calc_cond) { for (octave_idx_type j = 0; j < nr; j++) { double atmp = 0.; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) atmp += std::abs (data (i)); if (atmp > anorm) anorm = atmp; } } Array<octave_idx_type> ipvt (dim_vector (nr, 1)); octave_idx_type *pipvt = ipvt.fortran_vec (); F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt, err)); if (err != 0) { rcond = 0.0; err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (); } else { if (calc_cond) { char job = '1'; Array<Complex> z (dim_vector (2 * nr, 1)); Complex *pz = z.fortran_vec (); Array<double> iz (dim_vector (nr, 1)); double *piz = iz.fortran_vec (); F77_XFCN (zgbcon, ZGBCON, (F77_CONST_CHAR_ARG2 (&job, 1), nc, n_lower, n_upper, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt, anorm, rcond, F77_DBLE_CMPLX_ARG (pz), piz, err F77_CHAR_ARG_LEN (1))); if (err != 0) err = -2; volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } else rcond = 1.; if (err == 0) { char job = 'N'; volatile octave_idx_type x_nz = b.nnz (); octave_idx_type b_nc = b.cols (); retval = SparseComplexMatrix (nr, b_nc, x_nz); retval.xcidx (0) = 0; volatile octave_idx_type ii = 0; OCTAVE_LOCAL_BUFFER (Complex, work, nr); for (volatile octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) work[i] = 0.; for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) work[b.ridx (i)] = b.data (i); F77_XFCN (zgbtrs, ZGBTRS, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, n_upper, 1, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt, F77_DBLE_CMPLX_ARG (work), b.rows (), err F77_CHAR_ARG_LEN (1))); // Count nonzeros in work vector and adjust // space in retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nr; i++) if (work[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) { // Resize the sparse matrix octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; retval.change_capacity (sz); x_nz = sz; } for (octave_idx_type i = 0; i < nr; i++) if (work[i] != 0.) { retval.xridx (ii) = i; retval.xdata (ii++) = work[i]; } retval.xcidx (j+1) = ii; } retval.maybe_compress (); } } } else if (typ != MatrixType::Banded_Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; } ComplexMatrix SparseComplexMatrix::bsolve (MatrixType &mattype, const ComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); err = 0; if (nr != nc || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || b.cols () == 0) retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); else { // Print spparms("spumoni") info if requested volatile int typ = mattype.type (); mattype.info (); if (typ == MatrixType::Banded_Hermitian) { octave_idx_type n_lower = mattype.nlower (); octave_idx_type ldm = n_lower + 1; ComplexMatrix m_band (ldm, nc); Complex *tmp_data = m_band.fortran_vec (); if (! mattype.is_dense ()) { octave_idx_type ii = 0; for (octave_idx_type j = 0; j < ldm; j++) for (octave_idx_type i = 0; i < nc; i++) tmp_data[ii++] = 0.; } for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { octave_idx_type ri = ridx (i); if (ri >= j) m_band(ri - j, j) = data (i); } // Calculate the norm of the matrix, for later use. double anorm; if (calc_cond) anorm = m_band.abs ().sum ().row (0).max (); char job = 'L'; F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, F77_DBLE_CMPLX_ARG (tmp_data), ldm, err F77_CHAR_ARG_LEN (1))); if (err != 0) { // Matrix is not positive definite!! Fall through to // unsymmetric banded solver. rcond = 0.0; mattype.mark_as_unsymmetric (); typ = MatrixType::Banded; err = 0; } else { if (calc_cond) { Array<Complex> z (dim_vector (2 * nr, 1)); Complex *pz = z.fortran_vec (); Array<double> iz (dim_vector (nr, 1)); double *piz = iz.fortran_vec (); F77_XFCN (zpbcon, ZPBCON, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, F77_DBLE_CMPLX_ARG (tmp_data), ldm, anorm, rcond, F77_DBLE_CMPLX_ARG (pz), piz, err F77_CHAR_ARG_LEN (1))); if (err != 0) err = -2; volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } else rcond = 1.0; if (err == 0) { octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); retval = ComplexMatrix (b); Complex *result = retval.fortran_vec (); F77_XFCN (zpbtrs, ZPBTRS, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, b_nc, F77_DBLE_CMPLX_ARG (tmp_data), ldm, F77_DBLE_CMPLX_ARG (result), b_nr, err F77_CHAR_ARG_LEN (1))); if (err != 0) { // FIXME: Probably should be a warning. (*current_liboctave_error_handler) ("SparseComplexMatrix::solve solve failed"); err = -1; } } } } if (typ == MatrixType::Banded) { // Create the storage for the banded form of the sparse matrix octave_idx_type n_upper = mattype.nupper (); octave_idx_type n_lower = mattype.nlower (); octave_idx_type ldm = n_upper + 2 * n_lower + 1; ComplexMatrix m_band (ldm, nc); Complex *tmp_data = m_band.fortran_vec (); if (! mattype.is_dense ()) { octave_idx_type ii = 0; for (octave_idx_type j = 0; j < ldm; j++) for (octave_idx_type i = 0; i < nc; i++) tmp_data[ii++] = 0.; } for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) m_band(ridx (i) - j + n_lower + n_upper, j) = data (i); // Calculate the norm of the matrix, for later use. double anorm = 0.0; if (calc_cond) { for (octave_idx_type j = 0; j < nr; j++) { double atmp = 0.; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) atmp += std::abs (data (i)); if (atmp > anorm) anorm = atmp; } } Array<octave_idx_type> ipvt (dim_vector (nr, 1)); octave_idx_type *pipvt = ipvt.fortran_vec (); F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt, err)); if (err != 0) { err = -2; rcond = 0.0; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (); } else { if (calc_cond) { char job = '1'; Array<Complex> z (dim_vector (2 * nr, 1)); Complex *pz = z.fortran_vec (); Array<double> iz (dim_vector (nr, 1)); double *piz = iz.fortran_vec (); F77_XFCN (zgbcon, ZGBCON, (F77_CONST_CHAR_ARG2 (&job, 1), nc, n_lower, n_upper, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt, anorm, rcond, F77_DBLE_CMPLX_ARG (pz), piz, err F77_CHAR_ARG_LEN (1))); if (err != 0) err = -2; volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } else rcond = 1.; if (err == 0) { char job = 'N'; octave_idx_type b_nc = b.cols (); retval = ComplexMatrix (b); Complex *result = retval.fortran_vec (); F77_XFCN (zgbtrs, ZGBTRS, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, n_upper, b_nc, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt, F77_DBLE_CMPLX_ARG (result), b.rows (), err F77_CHAR_ARG_LEN (1))); } } } else if (typ != MatrixType::Banded_Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; } SparseComplexMatrix SparseComplexMatrix::bsolve (MatrixType &mattype, const SparseComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); err = 0; if (nr != nc || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || b.cols () == 0) retval = SparseComplexMatrix (nc, b.cols ()); else { // Print spparms("spumoni") info if requested volatile int typ = mattype.type (); mattype.info (); if (typ == MatrixType::Banded_Hermitian) { octave_idx_type n_lower = mattype.nlower (); octave_idx_type ldm = n_lower + 1; ComplexMatrix m_band (ldm, nc); Complex *tmp_data = m_band.fortran_vec (); if (! mattype.is_dense ()) { octave_idx_type ii = 0; for (octave_idx_type j = 0; j < ldm; j++) for (octave_idx_type i = 0; i < nc; i++) tmp_data[ii++] = 0.; } for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) { octave_idx_type ri = ridx (i); if (ri >= j) m_band(ri - j, j) = data (i); } // Calculate the norm of the matrix, for later use. double anorm; if (calc_cond) anorm = m_band.abs ().sum ().row (0).max (); char job = 'L'; F77_XFCN (zpbtrf, ZPBTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, F77_DBLE_CMPLX_ARG (tmp_data), ldm, err F77_CHAR_ARG_LEN (1))); if (err != 0) { // Matrix is not positive definite!! Fall through to // unsymmetric banded solver. mattype.mark_as_unsymmetric (); typ = MatrixType::Banded; rcond = 0.0; err = 0; } else { if (calc_cond) { Array<Complex> z (dim_vector (2 * nr, 1)); Complex *pz = z.fortran_vec (); Array<double> iz (dim_vector (nr, 1)); double *piz = iz.fortran_vec (); F77_XFCN (zpbcon, ZPBCON, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, F77_DBLE_CMPLX_ARG (tmp_data), ldm, anorm, rcond, F77_DBLE_CMPLX_ARG (pz), piz, err F77_CHAR_ARG_LEN (1))); if (err != 0) err = -2; volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } else rcond = 1.0; if (err == 0) { octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); // Take a first guess that the number of nonzero terms // will be as many as in b volatile octave_idx_type x_nz = b.nnz (); volatile octave_idx_type ii = 0; retval = SparseComplexMatrix (b_nr, b_nc, x_nz); retval.xcidx (0) = 0; for (volatile octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < b_nr; i++) Bx[i] = b(i,j); F77_XFCN (zpbtrs, ZPBTRS, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, 1, F77_DBLE_CMPLX_ARG (tmp_data), ldm, F77_DBLE_CMPLX_ARG (Bx), b_nr, err F77_CHAR_ARG_LEN (1))); if (err != 0) { // FIXME: Probably should be a warning. (*current_liboctave_error_handler) ("SparseMatrix::solve solve failed"); err = -1; break; } // Count nonzeros in work vector and adjust // space in retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nr; i++) if (Bx[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) { // Resize the sparse matrix octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; retval.change_capacity (sz); x_nz = sz; } for (octave_idx_type i = 0; i < nr; i++) if (Bx[i] != 0.) { retval.xridx (ii) = i; retval.xdata (ii++) = Bx[i]; } retval.xcidx (j+1) = ii; } retval.maybe_compress (); } } } if (typ == MatrixType::Banded) { // Create the storage for the banded form of the sparse matrix octave_idx_type n_upper = mattype.nupper (); octave_idx_type n_lower = mattype.nlower (); octave_idx_type ldm = n_upper + 2 * n_lower + 1; ComplexMatrix m_band (ldm, nc); Complex *tmp_data = m_band.fortran_vec (); if (! mattype.is_dense ()) { octave_idx_type ii = 0; for (octave_idx_type j = 0; j < ldm; j++) for (octave_idx_type i = 0; i < nc; i++) tmp_data[ii++] = 0.; } for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) m_band(ridx (i) - j + n_lower + n_upper, j) = data (i); // Calculate the norm of the matrix, for later use. double anorm = 0.0; if (calc_cond) { for (octave_idx_type j = 0; j < nr; j++) { double atmp = 0.; for (octave_idx_type i = cidx (j); i < cidx (j+1); i++) atmp += std::abs (data (i)); if (atmp > anorm) anorm = atmp; } } Array<octave_idx_type> ipvt (dim_vector (nr, 1)); octave_idx_type *pipvt = ipvt.fortran_vec (); F77_XFCN (zgbtrf, ZGBTRF, (nr, nr, n_lower, n_upper, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt, err)); if (err != 0) { err = -2; rcond = 0.0; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (); } else { if (calc_cond) { char job = '1'; Array<Complex> z (dim_vector (2 * nr, 1)); Complex *pz = z.fortran_vec (); Array<double> iz (dim_vector (nr, 1)); double *piz = iz.fortran_vec (); F77_XFCN (zgbcon, ZGBCON, (F77_CONST_CHAR_ARG2 (&job, 1), nc, n_lower, n_upper, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt, anorm, rcond, F77_DBLE_CMPLX_ARG (pz), piz, err F77_CHAR_ARG_LEN (1))); if (err != 0) err = -2; volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); } } else rcond = 1.; if (err == 0) { char job = 'N'; volatile octave_idx_type x_nz = b.nnz (); octave_idx_type b_nc = b.cols (); retval = SparseComplexMatrix (nr, b_nc, x_nz); retval.xcidx (0) = 0; volatile octave_idx_type ii = 0; OCTAVE_LOCAL_BUFFER (Complex, Bx, nr); for (volatile octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < nr; i++) Bx[i] = 0.; for (octave_idx_type i = b.cidx (j); i < b.cidx (j+1); i++) Bx[b.ridx (i)] = b.data (i); F77_XFCN (zgbtrs, ZGBTRS, (F77_CONST_CHAR_ARG2 (&job, 1), nr, n_lower, n_upper, 1, F77_DBLE_CMPLX_ARG (tmp_data), ldm, pipvt, F77_DBLE_CMPLX_ARG (Bx), b.rows (), err F77_CHAR_ARG_LEN (1))); // Count nonzeros in work vector and adjust // space in retval if needed octave_idx_type new_nnz = 0; for (octave_idx_type i = 0; i < nr; i++) if (Bx[i] != 0.) new_nnz++; if (ii + new_nnz > x_nz) { // Resize the sparse matrix octave_idx_type sz = new_nnz * (b_nc - j) + x_nz; retval.change_capacity (sz); x_nz = sz; } for (octave_idx_type i = 0; i < nr; i++) if (Bx[i] != 0.) { retval.xridx (ii) = i; retval.xdata (ii++) = Bx[i]; } retval.xcidx (j+1) = ii; } retval.maybe_compress (); } } } else if (typ != MatrixType::Banded_Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; } void * SparseComplexMatrix::factorize (octave_idx_type& err, double &rcond, Matrix &Control, Matrix &Info, solve_singularity_handler sing_handler, bool calc_cond) const { // The return values void *Numeric = 0; err = 0; #if defined (HAVE_UMFPACK) // Setup the control parameters Control = Matrix (UMFPACK_CONTROL, 1); double *control = Control.fortran_vec (); UMFPACK_ZNAME (defaults) (control); double tmp = octave_sparse_params::get_key ("spumoni"); if (! octave::math::isnan (tmp)) Control (UMFPACK_PRL) = tmp; tmp = octave_sparse_params::get_key ("piv_tol"); if (! octave::math::isnan (tmp)) { Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; Control (UMFPACK_PIVOT_TOLERANCE) = tmp; } // Set whether we are allowed to modify Q or not tmp = octave_sparse_params::get_key ("autoamd"); if (! octave::math::isnan (tmp)) Control (UMFPACK_FIXQ) = tmp; UMFPACK_ZNAME (report_control) (control); const octave_idx_type *Ap = cidx (); const octave_idx_type *Ai = ridx (); const Complex *Ax = data (); octave_idx_type nr = rows (); octave_idx_type nc = cols (); UMFPACK_ZNAME (report_matrix) (nr, nc, Ap, Ai, reinterpret_cast<const double *> (Ax), 0, 1, control); void *Symbolic; Info = Matrix (1, UMFPACK_INFO); double *info = Info.fortran_vec (); int status = UMFPACK_ZNAME (qsymbolic) (nr, nc, Ap, Ai, reinterpret_cast<const double *> (Ax), 0, 0, &Symbolic, control, info); if (status < 0) { UMFPACK_ZNAME (report_status) (control, status); UMFPACK_ZNAME (report_info) (control, info); UMFPACK_ZNAME (free_symbolic) (&Symbolic); // FIXME: Should this be a warning? (*current_liboctave_error_handler) ("SparseComplexMatrix::solve symbolic factorization failed"); err = -1; } else { UMFPACK_ZNAME (report_symbolic) (Symbolic, control); status = UMFPACK_ZNAME (numeric) (Ap, Ai, reinterpret_cast<const double *> (Ax), 0, Symbolic, &Numeric, control, info); UMFPACK_ZNAME (free_symbolic) (&Symbolic); if (calc_cond) rcond = Info (UMFPACK_RCOND); else rcond = 1.; volatile double rcond_plus_one = rcond + 1.0; if (status == UMFPACK_WARNING_singular_matrix || rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { UMFPACK_ZNAME (report_numeric) (Numeric, control); err = -2; if (sing_handler) sing_handler (rcond); else octave::warn_singular_matrix (rcond); } else if (status < 0) { UMFPACK_ZNAME (report_status) (control, status); UMFPACK_ZNAME (report_info) (control, info); // FIXME: Should this be a warning? (*current_liboctave_error_handler) ("SparseComplexMatrix::solve numeric factorization failed"); err = -1; } else { UMFPACK_ZNAME (report_numeric) (Numeric, control); } } if (err != 0) UMFPACK_ZNAME (free_numeric) (&Numeric); #else octave_unused_parameter (rcond); octave_unused_parameter (Control); octave_unused_parameter (Info); octave_unused_parameter (sing_handler); octave_unused_parameter (calc_cond); (*current_liboctave_error_handler) ("support for UMFPACK was unavailable or disabled when liboctave was built"); #endif return Numeric; } ComplexMatrix SparseComplexMatrix::fsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); err = 0; if (nr != nc || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || b.cols () == 0) retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); else { // Print spparms("spumoni") info if requested volatile int typ = mattype.type (); mattype.info (); if (typ == MatrixType::Hermitian) { #if defined (HAVE_CHOLMOD) cholmod_common Common; cholmod_common *cm = &Common; // Setup initial parameters CHOLMOD_NAME(start) (cm); cm->prefer_zomplex = false; double spu = octave_sparse_params::get_key ("spumoni"); if (spu == 0.) { cm->print = -1; SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0); } else { cm->print = static_cast<int> (spu) + 2; SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint); } cm->error_handler = &SparseCholError; SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex); SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot); cm->final_ll = true; cholmod_sparse Astore; cholmod_sparse *A = &Astore; double dummy; A->nrow = nr; A->ncol = nc; A->p = cidx (); A->i = ridx (); A->nzmax = nnz (); A->packed = true; A->sorted = true; A->nz = 0; #if defined (OCTAVE_ENABLE_64) A->itype = CHOLMOD_LONG; #else A->itype = CHOLMOD_INT; #endif A->dtype = CHOLMOD_DOUBLE; A->stype = 1; A->xtype = CHOLMOD_COMPLEX; if (nr < 1) A->x = &dummy; else A->x = data (); cholmod_dense Bstore; cholmod_dense *B = &Bstore; B->nrow = b.rows (); B->ncol = b.cols (); B->d = B->nrow; B->nzmax = B->nrow * B->ncol; B->dtype = CHOLMOD_DOUBLE; B->xtype = CHOLMOD_REAL; if (nc < 1 || b.cols () < 1) B->x = &dummy; else // We won't alter it, honest :-) B->x = const_cast<double *>(b.fortran_vec ()); cholmod_factor *L; BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; L = CHOLMOD_NAME(analyze) (A, cm); CHOLMOD_NAME(factorize) (A, L, cm); if (calc_cond) rcond = CHOLMOD_NAME(rcond)(L, cm); else rcond = 1.; END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; if (rcond == 0.0) { // Either its indefinite or singular. Try UMFPACK mattype.mark_as_unsymmetric (); typ = MatrixType::Full; } else { volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); return retval; } cholmod_dense *X; BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm); END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; retval.resize (b.rows (), b.cols ()); for (octave_idx_type j = 0; j < b.cols (); j++) { octave_idx_type jr = j * b.rows (); for (octave_idx_type i = 0; i < b.rows (); i++) retval.xelem (i,j) = static_cast<Complex *>(X->x)[jr + i]; } BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; CHOLMOD_NAME(free_dense) (&X, cm); CHOLMOD_NAME(free_factor) (&L, cm); CHOLMOD_NAME(finish) (cm); static char tmp[] = " "; CHOLMOD_NAME(print_common) (tmp, cm); END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; } #else (*current_liboctave_warning_with_id_handler) ("Octave:missing-dependency", "support for CHOLMOD was unavailable or disabled " "when liboctave was built"); mattype.mark_as_unsymmetric (); typ = MatrixType::Full; #endif } if (typ == MatrixType::Full) { #if defined (HAVE_UMFPACK) Matrix Control, Info; void *Numeric = factorize (err, rcond, Control, Info, sing_handler, calc_cond); if (err == 0) { octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); int status = 0; double *control = Control.fortran_vec (); double *info = Info.fortran_vec (); const octave_idx_type *Ap = cidx (); const octave_idx_type *Ai = ridx (); const Complex *Ax = data (); #if defined (UMFPACK_SEPARATE_SPLIT) const double *Bx = b.fortran_vec (); OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); for (octave_idx_type i = 0; i < b_nr; i++) Bz[i] = 0.; #else OCTAVE_LOCAL_BUFFER (Complex, Bz, b_nr); #endif retval.resize (b_nr, b_nc); Complex *Xx = retval.fortran_vec (); for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr) { #if defined (UMFPACK_SEPARATE_SPLIT) status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai, reinterpret_cast<const double *> (Ax), 0, reinterpret_cast<double *> (&Xx[iidx]), 0, &Bx[iidx], Bz, Numeric, control, info); #else for (octave_idx_type i = 0; i < b_nr; i++) Bz[i] = b.elem (i, j); status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai, reinterpret_cast<const double *> (Ax), 0, reinterpret_cast<double *> (&Xx[iidx]), 0, reinterpret_cast<const double *> (Bz), 0, Numeric, control, info); #endif if (status < 0) { UMFPACK_ZNAME (report_status) (control, status); // FIXME: Should this be a warning? (*current_liboctave_error_handler) ("SparseComplexMatrix::solve solve failed"); err = -1; break; } } UMFPACK_ZNAME (report_info) (control, info); UMFPACK_ZNAME (free_numeric) (&Numeric); } else mattype.mark_as_rectangular (); #else octave_unused_parameter (rcond); octave_unused_parameter (sing_handler); octave_unused_parameter (calc_cond); (*current_liboctave_error_handler) ("support for UMFPACK was unavailable or disabled " "when liboctave was built"); #endif } else if (typ != MatrixType::Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; } SparseComplexMatrix SparseComplexMatrix::fsolve (MatrixType &mattype, const SparseMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); err = 0; if (nr != nc || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || b.cols () == 0) retval = SparseComplexMatrix (nc, b.cols ()); else { // Print spparms("spumoni") info if requested volatile int typ = mattype.type (); mattype.info (); if (typ == MatrixType::Hermitian) { #if defined (HAVE_CHOLMOD) cholmod_common Common; cholmod_common *cm = &Common; // Setup initial parameters CHOLMOD_NAME(start) (cm); cm->prefer_zomplex = false; double spu = octave_sparse_params::get_key ("spumoni"); if (spu == 0.) { cm->print = -1; SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0); } else { cm->print = static_cast<int> (spu) + 2; SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint); } cm->error_handler = &SparseCholError; SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex); SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot); cm->final_ll = true; cholmod_sparse Astore; cholmod_sparse *A = &Astore; double dummy; A->nrow = nr; A->ncol = nc; A->p = cidx (); A->i = ridx (); A->nzmax = nnz (); A->packed = true; A->sorted = true; A->nz = 0; #if defined (OCTAVE_ENABLE_64) A->itype = CHOLMOD_LONG; #else A->itype = CHOLMOD_INT; #endif A->dtype = CHOLMOD_DOUBLE; A->stype = 1; A->xtype = CHOLMOD_COMPLEX; if (nr < 1) A->x = &dummy; else A->x = data (); cholmod_sparse Bstore; cholmod_sparse *B = &Bstore; B->nrow = b.rows (); B->ncol = b.cols (); B->p = b.cidx (); B->i = b.ridx (); B->nzmax = b.nnz (); B->packed = true; B->sorted = true; B->nz = 0; #if defined (OCTAVE_ENABLE_64) B->itype = CHOLMOD_LONG; #else B->itype = CHOLMOD_INT; #endif B->dtype = CHOLMOD_DOUBLE; B->stype = 0; B->xtype = CHOLMOD_REAL; if (b.rows () < 1 || b.cols () < 1) B->x = &dummy; else B->x = b.data (); cholmod_factor *L; BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; L = CHOLMOD_NAME(analyze) (A, cm); CHOLMOD_NAME(factorize) (A, L, cm); if (calc_cond) rcond = CHOLMOD_NAME(rcond)(L, cm); else rcond = 1.; END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; if (rcond == 0.0) { // Either its indefinite or singular. Try UMFPACK mattype.mark_as_unsymmetric (); typ = MatrixType::Full; } else { volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); return retval; } cholmod_sparse *X; BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm); END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; retval = SparseComplexMatrix (static_cast<octave_idx_type>(X->nrow), static_cast<octave_idx_type>(X->ncol), static_cast<octave_idx_type>(X->nzmax)); for (octave_idx_type j = 0; j <= static_cast<octave_idx_type>(X->ncol); j++) retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j]; for (octave_idx_type j = 0; j < static_cast<octave_idx_type>(X->nzmax); j++) { retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j]; retval.xdata (j) = static_cast<Complex *>(X->x)[j]; } BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; CHOLMOD_NAME(free_sparse) (&X, cm); CHOLMOD_NAME(free_factor) (&L, cm); CHOLMOD_NAME(finish) (cm); static char tmp[] = " "; CHOLMOD_NAME(print_common) (tmp, cm); END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; } #else (*current_liboctave_warning_with_id_handler) ("Octave:missing-dependency", "support for CHOLMOD was unavailable or disabled " "when liboctave was built"); mattype.mark_as_unsymmetric (); typ = MatrixType::Full; #endif } if (typ == MatrixType::Full) { #if defined (HAVE_UMFPACK) Matrix Control, Info; void *Numeric = factorize (err, rcond, Control, Info, sing_handler, calc_cond); if (err == 0) { octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); int status = 0; double *control = Control.fortran_vec (); double *info = Info.fortran_vec (); const octave_idx_type *Ap = cidx (); const octave_idx_type *Ai = ridx (); const Complex *Ax = data (); #if defined (UMFPACK_SEPARATE_SPLIT) OCTAVE_LOCAL_BUFFER (double, Bx, b_nr); OCTAVE_LOCAL_BUFFER (double, Bz, b_nr); for (octave_idx_type i = 0; i < b_nr; i++) Bz[i] = 0.; #else OCTAVE_LOCAL_BUFFER (Complex, Bz, b_nr); #endif // Take a first guess that the number of nonzero terms // will be as many as in b octave_idx_type x_nz = b.nnz (); octave_idx_type ii = 0; retval = SparseComplexMatrix (b_nr, b_nc, x_nz); OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr); retval.xcidx (0) = 0; for (octave_idx_type j = 0; j < b_nc; j++) { #if defined (UMFPACK_SEPARATE_SPLIT) for (octave_idx_type i = 0; i < b_nr; i++) Bx[i] = b.elem (i, j); status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai, reinterpret_cast<const double *> (Ax), 0, reinterpret_cast<double *> (Xx), 0, Bx, Bz, Numeric, control, info); #else for (octave_idx_type i = 0; i < b_nr; i++) Bz[i] = b.elem (i, j); status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai, reinterpret_cast<const double *> (Ax), 0, reinterpret_cast<double *> (Xx), 0, reinterpret_cast<double *> (Bz), 0, Numeric, control, info); #endif if (status < 0) { UMFPACK_ZNAME (report_status) (control, status); // FIXME: Should this be a warning? (*current_liboctave_error_handler) ("SparseComplexMatrix::solve solve failed"); err = -1; break; } for (octave_idx_type i = 0; i < b_nr; i++) { Complex tmp = Xx[i]; if (tmp != 0.0) { if (ii == x_nz) { // Resize the sparse matrix octave_idx_type sz = x_nz * (b_nc - j) / b_nc; sz = (sz > 10 ? sz : 10) + x_nz; retval.change_capacity (sz); x_nz = sz; } retval.xdata (ii) = tmp; retval.xridx (ii++) = i; } } retval.xcidx (j+1) = ii; } retval.maybe_compress (); UMFPACK_ZNAME (report_info) (control, info); UMFPACK_ZNAME (free_numeric) (&Numeric); } else mattype.mark_as_rectangular (); #else octave_unused_parameter (rcond); octave_unused_parameter (sing_handler); octave_unused_parameter (calc_cond); (*current_liboctave_error_handler) ("support for UMFPACK was unavailable or disabled " "when liboctave was built"); #endif } else if (typ != MatrixType::Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; } ComplexMatrix SparseComplexMatrix::fsolve (MatrixType &mattype, const ComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { ComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); err = 0; if (nr != nc || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || b.cols () == 0) retval = ComplexMatrix (nc, b.cols (), Complex (0.0, 0.0)); else { // Print spparms("spumoni") info if requested volatile int typ = mattype.type (); mattype.info (); if (typ == MatrixType::Hermitian) { #if defined (HAVE_CHOLMOD) cholmod_common Common; cholmod_common *cm = &Common; // Setup initial parameters CHOLMOD_NAME(start) (cm); cm->prefer_zomplex = false; double spu = octave_sparse_params::get_key ("spumoni"); if (spu == 0.) { cm->print = -1; SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0); } else { cm->print = static_cast<int> (spu) + 2; SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint); } cm->error_handler = &SparseCholError; SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex); SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot); cm->final_ll = true; cholmod_sparse Astore; cholmod_sparse *A = &Astore; double dummy; A->nrow = nr; A->ncol = nc; A->p = cidx (); A->i = ridx (); A->nzmax = nnz (); A->packed = true; A->sorted = true; A->nz = 0; #if defined (OCTAVE_ENABLE_64) A->itype = CHOLMOD_LONG; #else A->itype = CHOLMOD_INT; #endif A->dtype = CHOLMOD_DOUBLE; A->stype = 1; A->xtype = CHOLMOD_COMPLEX; if (nr < 1) A->x = &dummy; else A->x = data (); cholmod_dense Bstore; cholmod_dense *B = &Bstore; B->nrow = b.rows (); B->ncol = b.cols (); B->d = B->nrow; B->nzmax = B->nrow * B->ncol; B->dtype = CHOLMOD_DOUBLE; B->xtype = CHOLMOD_COMPLEX; if (nc < 1 || b.cols () < 1) B->x = &dummy; else // We won't alter it, honest :-) B->x = const_cast<Complex *>(b.fortran_vec ()); cholmod_factor *L; BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; L = CHOLMOD_NAME(analyze) (A, cm); CHOLMOD_NAME(factorize) (A, L, cm); if (calc_cond) rcond = CHOLMOD_NAME(rcond)(L, cm); else rcond = 1.; END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; if (rcond == 0.0) { // Either its indefinite or singular. Try UMFPACK mattype.mark_as_unsymmetric (); typ = MatrixType::Full; } else { volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); return retval; } cholmod_dense *X; BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; X = CHOLMOD_NAME(solve) (CHOLMOD_A, L, B, cm); END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; retval.resize (b.rows (), b.cols ()); for (octave_idx_type j = 0; j < b.cols (); j++) { octave_idx_type jr = j * b.rows (); for (octave_idx_type i = 0; i < b.rows (); i++) retval.xelem (i,j) = static_cast<Complex *>(X->x)[jr + i]; } BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; CHOLMOD_NAME(free_dense) (&X, cm); CHOLMOD_NAME(free_factor) (&L, cm); CHOLMOD_NAME(finish) (cm); static char tmp[] = " "; CHOLMOD_NAME(print_common) (tmp, cm); END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; } #else (*current_liboctave_warning_with_id_handler) ("Octave:missing-dependency", "support for CHOLMOD was unavailable or disabled " "when liboctave was built"); mattype.mark_as_unsymmetric (); typ = MatrixType::Full; #endif } if (typ == MatrixType::Full) { #if defined (HAVE_UMFPACK) Matrix Control, Info; void *Numeric = factorize (err, rcond, Control, Info, sing_handler, calc_cond); if (err == 0) { octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); int status = 0; double *control = Control.fortran_vec (); double *info = Info.fortran_vec (); const octave_idx_type *Ap = cidx (); const octave_idx_type *Ai = ridx (); const Complex *Ax = data (); const Complex *Bx = b.fortran_vec (); retval.resize (b_nr, b_nc); Complex *Xx = retval.fortran_vec (); for (octave_idx_type j = 0, iidx = 0; j < b_nc; j++, iidx += b_nr) { status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai, reinterpret_cast<const double *> (Ax), 0, reinterpret_cast<double *> (&Xx[iidx]), 0, reinterpret_cast<const double *> (&Bx[iidx]), 0, Numeric, control, info); if (status < 0) { UMFPACK_ZNAME (report_status) (control, status); // FIXME: Should this be a warning? (*current_liboctave_error_handler) ("SparseComplexMatrix::solve solve failed"); err = -1; break; } } UMFPACK_ZNAME (report_info) (control, info); UMFPACK_ZNAME (free_numeric) (&Numeric); } else mattype.mark_as_rectangular (); #else octave_unused_parameter (rcond); octave_unused_parameter (sing_handler); octave_unused_parameter (calc_cond); (*current_liboctave_error_handler) ("support for UMFPACK was unavailable or disabled " "when liboctave was built"); #endif } else if (typ != MatrixType::Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; } SparseComplexMatrix SparseComplexMatrix::fsolve (MatrixType &mattype, const SparseComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool calc_cond) const { SparseComplexMatrix retval; octave_idx_type nr = rows (); octave_idx_type nc = cols (); err = 0; if (nr != nc || nr != b.rows ()) (*current_liboctave_error_handler) ("matrix dimension mismatch solution of linear equations"); if (nr == 0 || b.cols () == 0) retval = SparseComplexMatrix (nc, b.cols ()); else { // Print spparms("spumoni") info if requested volatile int typ = mattype.type (); mattype.info (); if (typ == MatrixType::Hermitian) { #if defined (HAVE_CHOLMOD) cholmod_common Common; cholmod_common *cm = &Common; // Setup initial parameters CHOLMOD_NAME(start) (cm); cm->prefer_zomplex = false; double spu = octave_sparse_params::get_key ("spumoni"); if (spu == 0.) { cm->print = -1; SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0); } else { cm->print = static_cast<int> (spu) + 2; SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint); } cm->error_handler = &SparseCholError; SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex); SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot); cm->final_ll = true; cholmod_sparse Astore; cholmod_sparse *A = &Astore; double dummy; A->nrow = nr; A->ncol = nc; A->p = cidx (); A->i = ridx (); A->nzmax = nnz (); A->packed = true; A->sorted = true; A->nz = 0; #if defined (OCTAVE_ENABLE_64) A->itype = CHOLMOD_LONG; #else A->itype = CHOLMOD_INT; #endif A->dtype = CHOLMOD_DOUBLE; A->stype = 1; A->xtype = CHOLMOD_COMPLEX; if (nr < 1) A->x = &dummy; else A->x = data (); cholmod_sparse Bstore; cholmod_sparse *B = &Bstore; B->nrow = b.rows (); B->ncol = b.cols (); B->p = b.cidx (); B->i = b.ridx (); B->nzmax = b.nnz (); B->packed = true; B->sorted = true; B->nz = 0; #if defined (OCTAVE_ENABLE_64) B->itype = CHOLMOD_LONG; #else B->itype = CHOLMOD_INT; #endif B->dtype = CHOLMOD_DOUBLE; B->stype = 0; B->xtype = CHOLMOD_COMPLEX; if (b.rows () < 1 || b.cols () < 1) B->x = &dummy; else B->x = b.data (); cholmod_factor *L; BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; L = CHOLMOD_NAME(analyze) (A, cm); CHOLMOD_NAME(factorize) (A, L, cm); if (calc_cond) rcond = CHOLMOD_NAME(rcond)(L, cm); else rcond = 1.; END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; if (rcond == 0.0) { // Either its indefinite or singular. Try UMFPACK mattype.mark_as_unsymmetric (); typ = MatrixType::Full; } else { volatile double rcond_plus_one = rcond + 1.0; if (rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) { sing_handler (rcond); mattype.mark_as_rectangular (); } else octave::warn_singular_matrix (rcond); return retval; } cholmod_sparse *X; BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; X = CHOLMOD_NAME(spsolve) (CHOLMOD_A, L, B, cm); END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; retval = SparseComplexMatrix (static_cast<octave_idx_type>(X->nrow), static_cast<octave_idx_type>(X->ncol), static_cast<octave_idx_type>(X->nzmax)); for (octave_idx_type j = 0; j <= static_cast<octave_idx_type>(X->ncol); j++) retval.xcidx (j) = static_cast<octave_idx_type *>(X->p)[j]; for (octave_idx_type j = 0; j < static_cast<octave_idx_type>(X->nzmax); j++) { retval.xridx (j) = static_cast<octave_idx_type *>(X->i)[j]; retval.xdata (j) = static_cast<Complex *>(X->x)[j]; } BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; CHOLMOD_NAME(free_sparse) (&X, cm); CHOLMOD_NAME(free_factor) (&L, cm); CHOLMOD_NAME(finish) (cm); static char tmp[] = " "; CHOLMOD_NAME(print_common) (tmp, cm); END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; } #else (*current_liboctave_warning_with_id_handler) ("Octave:missing-dependency", "support for CHOLMOD was unavailable or disabled " "when liboctave was built"); mattype.mark_as_unsymmetric (); typ = MatrixType::Full; #endif } if (typ == MatrixType::Full) { #if defined (HAVE_UMFPACK) Matrix Control, Info; void *Numeric = factorize (err, rcond, Control, Info, sing_handler, calc_cond); if (err == 0) { octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); int status = 0; double *control = Control.fortran_vec (); double *info = Info.fortran_vec (); const octave_idx_type *Ap = cidx (); const octave_idx_type *Ai = ridx (); const Complex *Ax = data (); OCTAVE_LOCAL_BUFFER (Complex, Bx, b_nr); // Take a first guess that the number of nonzero terms // will be as many as in b octave_idx_type x_nz = b.nnz (); octave_idx_type ii = 0; retval = SparseComplexMatrix (b_nr, b_nc, x_nz); OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr); retval.xcidx (0) = 0; for (octave_idx_type j = 0; j < b_nc; j++) { for (octave_idx_type i = 0; i < b_nr; i++) Bx[i] = b(i,j); status = UMFPACK_ZNAME (solve) (UMFPACK_A, Ap, Ai, reinterpret_cast<const double *> (Ax), 0, reinterpret_cast<double *> (Xx), 0, reinterpret_cast<double *> (Bx), 0, Numeric, control, info); if (status < 0) { UMFPACK_ZNAME (report_status) (control, status); // FIXME: Should this be a warning? (*current_liboctave_error_handler) ("SparseComplexMatrix::solve solve failed"); err = -1; break; } for (octave_idx_type i = 0; i < b_nr; i++) { Complex tmp = Xx[i]; if (tmp != 0.0) { if (ii == x_nz) { // Resize the sparse matrix octave_idx_type sz = x_nz * (b_nc - j) / b_nc; sz = (sz > 10 ? sz : 10) + x_nz; retval.change_capacity (sz); x_nz = sz; } retval.xdata (ii) = tmp; retval.xridx (ii++) = i; } } retval.xcidx (j+1) = ii; } retval.maybe_compress (); rcond = Info (UMFPACK_RCOND); volatile double rcond_plus_one = rcond + 1.0; if (status == UMFPACK_WARNING_singular_matrix || rcond_plus_one == 1.0 || octave::math::isnan (rcond)) { err = -2; if (sing_handler) sing_handler (rcond); else octave::warn_singular_matrix (rcond); } UMFPACK_ZNAME (report_info) (control, info); UMFPACK_ZNAME (free_numeric) (&Numeric); } else mattype.mark_as_rectangular (); #else octave_unused_parameter (rcond); octave_unused_parameter (sing_handler); octave_unused_parameter (calc_cond); (*current_liboctave_error_handler) ("support for UMFPACK was unavailable or disabled " "when liboctave was built"); #endif } else if (typ != MatrixType::Hermitian) (*current_liboctave_error_handler) ("incorrect matrix type"); } return retval; } ComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const Matrix& b) const { octave_idx_type info; double rcond; return solve (mattype, b, info, rcond, 0); } ComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const Matrix& b, octave_idx_type& info) const { double rcond; return solve (mattype, b, info, rcond, 0); } ComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const Matrix& b, octave_idx_type& info, double& rcond) const { return solve (mattype, b, info, rcond, 0); } ComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const Matrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool singular_fallback) const { ComplexMatrix retval; int typ = mattype.type (false); if (typ == MatrixType::Unknown) typ = mattype.type (*this); if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal) retval = dsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) retval = utsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) retval = ltsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian) retval = bsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Tridiagonal || typ == MatrixType::Tridiagonal_Hermitian) retval = trisolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) retval = fsolve (mattype, b, err, rcond, sing_handler, true); else if (typ != MatrixType::Rectangular) (*current_liboctave_error_handler) ("unknown matrix type"); if (singular_fallback && mattype.type (false) == MatrixType::Rectangular) { rcond = 1.; #if defined (USE_QRSOLVE) retval = qrsolve (*this, b, err); #else retval = dmsolve<ComplexMatrix, SparseComplexMatrix, Matrix> (*this, b, err); #endif } return retval; } SparseComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b) const { octave_idx_type info; double rcond; return solve (mattype, b, info, rcond, 0); } SparseComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b, octave_idx_type& info) const { double rcond; return solve (mattype, b, info, rcond, 0); } SparseComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b, octave_idx_type& info, double& rcond) const { return solve (mattype, b, info, rcond, 0); } SparseComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const SparseMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool singular_fallback) const { SparseComplexMatrix retval; int typ = mattype.type (false); if (typ == MatrixType::Unknown) typ = mattype.type (*this); if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal) retval = dsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) retval = utsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) retval = ltsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian) retval = bsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Tridiagonal || typ == MatrixType::Tridiagonal_Hermitian) retval = trisolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) retval = fsolve (mattype, b, err, rcond, sing_handler, true); else if (typ != MatrixType::Rectangular) (*current_liboctave_error_handler) ("unknown matrix type"); if (singular_fallback && mattype.type (false) == MatrixType::Rectangular) { rcond = 1.; #if defined (USE_QRSOLVE) retval = qrsolve (*this, b, err); #else retval = dmsolve<SparseComplexMatrix, SparseComplexMatrix, SparseMatrix> (*this, b, err); #endif } return retval; } ComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b) const { octave_idx_type info; double rcond; return solve (mattype, b, info, rcond, 0); } ComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b, octave_idx_type& info) const { double rcond; return solve (mattype, b, info, rcond, 0); } ComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b, octave_idx_type& info, double& rcond) const { return solve (mattype, b, info, rcond, 0); } ComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool singular_fallback) const { ComplexMatrix retval; int typ = mattype.type (false); if (typ == MatrixType::Unknown) typ = mattype.type (*this); if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal) retval = dsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) retval = utsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) retval = ltsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian) retval = bsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Tridiagonal || typ == MatrixType::Tridiagonal_Hermitian) retval = trisolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) retval = fsolve (mattype, b, err, rcond, sing_handler, true); else if (typ != MatrixType::Rectangular) (*current_liboctave_error_handler) ("unknown matrix type"); if (singular_fallback && mattype.type (false) == MatrixType::Rectangular) { rcond = 1.; #if defined (USE_QRSOLVE) retval = qrsolve (*this, b, err); #else retval = dmsolve<ComplexMatrix, SparseComplexMatrix, ComplexMatrix> (*this, b, err); #endif } return retval; } SparseComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b) const { octave_idx_type info; double rcond; return solve (mattype, b, info, rcond, 0); } SparseComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b, octave_idx_type& info) const { double rcond; return solve (mattype, b, info, rcond, 0); } SparseComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b, octave_idx_type& info, double& rcond) const { return solve (mattype, b, info, rcond, 0); } SparseComplexMatrix SparseComplexMatrix::solve (MatrixType &mattype, const SparseComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler, bool singular_fallback) const { SparseComplexMatrix retval; int typ = mattype.type (false); if (typ == MatrixType::Unknown) typ = mattype.type (*this); if (typ == MatrixType::Diagonal || typ == MatrixType::Permuted_Diagonal) retval = dsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper) retval = utsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower) retval = ltsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Banded || typ == MatrixType::Banded_Hermitian) retval = bsolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Tridiagonal || typ == MatrixType::Tridiagonal_Hermitian) retval = trisolve (mattype, b, err, rcond, sing_handler, false); else if (typ == MatrixType::Full || typ == MatrixType::Hermitian) retval = fsolve (mattype, b, err, rcond, sing_handler, true); else if (typ != MatrixType::Rectangular) (*current_liboctave_error_handler) ("unknown matrix type"); if (singular_fallback && mattype.type (false) == MatrixType::Rectangular) { rcond = 1.; #if defined (USE_QRSOLVE) retval = qrsolve (*this, b, err); #else retval = dmsolve<SparseComplexMatrix, SparseComplexMatrix, SparseComplexMatrix> (*this, b, err); #endif } return retval; } ComplexColumnVector SparseComplexMatrix::solve (MatrixType &mattype, const ColumnVector& b) const { octave_idx_type info; double rcond; return solve (mattype, b, info, rcond); } ComplexColumnVector SparseComplexMatrix::solve (MatrixType &mattype, const ColumnVector& b, octave_idx_type& info) const { double rcond; return solve (mattype, b, info, rcond); } ComplexColumnVector SparseComplexMatrix::solve (MatrixType &mattype, const ColumnVector& b, octave_idx_type& info, double& rcond) const { return solve (mattype, b, info, rcond, 0); } ComplexColumnVector SparseComplexMatrix::solve (MatrixType &mattype, const ColumnVector& b, octave_idx_type& info, double& rcond, solve_singularity_handler sing_handler) const { Matrix tmp (b); return solve (mattype, tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0)); } ComplexColumnVector SparseComplexMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b) const { octave_idx_type info; double rcond; return solve (mattype, b, info, rcond, 0); } ComplexColumnVector SparseComplexMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b, octave_idx_type& info) const { double rcond; return solve (mattype, b, info, rcond, 0); } ComplexColumnVector SparseComplexMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b, octave_idx_type& info, double& rcond) const { return solve (mattype, b, info, rcond, 0); } ComplexColumnVector SparseComplexMatrix::solve (MatrixType &mattype, const ComplexColumnVector& b, octave_idx_type& info, double& rcond, solve_singularity_handler sing_handler) const { ComplexMatrix tmp (b); return solve (mattype, tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0)); } ComplexMatrix SparseComplexMatrix::solve (const Matrix& b) const { octave_idx_type info; double rcond; return solve (b, info, rcond, 0); } ComplexMatrix SparseComplexMatrix::solve (const Matrix& b, octave_idx_type& info) const { double rcond; return solve (b, info, rcond, 0); } ComplexMatrix SparseComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcond) const { return solve (b, info, rcond, 0); } ComplexMatrix SparseComplexMatrix::solve (const Matrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler) const { MatrixType mattype (*this); return solve (mattype, b, err, rcond, sing_handler); } SparseComplexMatrix SparseComplexMatrix::solve (const SparseMatrix& b) const { octave_idx_type info; double rcond; return solve (b, info, rcond, 0); } SparseComplexMatrix SparseComplexMatrix::solve (const SparseMatrix& b, octave_idx_type& info) const { double rcond; return solve (b, info, rcond, 0); } SparseComplexMatrix SparseComplexMatrix::solve (const SparseMatrix& b, octave_idx_type& info, double& rcond) const { return solve (b, info, rcond, 0); } SparseComplexMatrix SparseComplexMatrix::solve (const SparseMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler) const { MatrixType mattype (*this); return solve (mattype, b, err, rcond, sing_handler); } ComplexMatrix SparseComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info) const { double rcond; return solve (b, info, rcond, 0); } ComplexMatrix SparseComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double& rcond) const { return solve (b, info, rcond, 0); } ComplexMatrix SparseComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler) const { MatrixType mattype (*this); return solve (mattype, b, err, rcond, sing_handler); } SparseComplexMatrix SparseComplexMatrix::solve (const SparseComplexMatrix& b) const { octave_idx_type info; double rcond; return solve (b, info, rcond, 0); } SparseComplexMatrix SparseComplexMatrix::solve (const SparseComplexMatrix& b, octave_idx_type& info) const { double rcond; return solve (b, info, rcond, 0); } SparseComplexMatrix SparseComplexMatrix::solve (const SparseComplexMatrix& b, octave_idx_type& info, double& rcond) const { return solve (b, info, rcond, 0); } SparseComplexMatrix SparseComplexMatrix::solve (const SparseComplexMatrix& b, octave_idx_type& err, double& rcond, solve_singularity_handler sing_handler) const { MatrixType mattype (*this); return solve (mattype, b, err, rcond, sing_handler); } ComplexColumnVector SparseComplexMatrix::solve (const ColumnVector& b) const { octave_idx_type info; double rcond; return solve (b, info, rcond); } ComplexColumnVector SparseComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info) const { double rcond; return solve (b, info, rcond); } ComplexColumnVector SparseComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcond) const { return solve (b, info, rcond, 0); } ComplexColumnVector SparseComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcond, solve_singularity_handler sing_handler) const { Matrix tmp (b); return solve (tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0)); } ComplexColumnVector SparseComplexMatrix::solve (const ComplexColumnVector& b) const { octave_idx_type info; double rcond; return solve (b, info, rcond, 0); } ComplexColumnVector SparseComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info) const { double rcond; return solve (b, info, rcond, 0); } ComplexColumnVector SparseComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info, double& rcond) const { return solve (b, info, rcond, 0); } ComplexColumnVector SparseComplexMatrix::solve (const ComplexColumnVector& b, octave_idx_type& info, double& rcond, solve_singularity_handler sing_handler) const { ComplexMatrix tmp (b); return solve (tmp, info, rcond, sing_handler).column (static_cast<octave_idx_type> (0)); } // unary operations SparseBoolMatrix SparseComplexMatrix::operator ! (void) const { if (any_element_is_nan ()) octave::err_nan_to_logical_conversion (); octave_idx_type nr = rows (); octave_idx_type nc = cols (); octave_idx_type nz1 = nnz (); octave_idx_type nz2 = nr*nc - nz1; SparseBoolMatrix r (nr, nc, nz2); octave_idx_type ii = 0; octave_idx_type jj = 0; r.cidx (0) = 0; for (octave_idx_type i = 0; i < nc; i++) { for (octave_idx_type j = 0; j < nr; j++) { if (jj < cidx (i+1) && ridx (jj) == j) jj++; else { r.data (ii) = true; r.ridx (ii++) = j; } } r.cidx (i+1) = ii; } return r; } SparseComplexMatrix SparseComplexMatrix::squeeze (void) const { return MSparse<Complex>::squeeze (); } SparseComplexMatrix SparseComplexMatrix::reshape (const dim_vector& new_dims) const { return MSparse<Complex>::reshape (new_dims); } SparseComplexMatrix SparseComplexMatrix::permute (const Array<octave_idx_type>& vec, bool inv) const { return MSparse<Complex>::permute (vec, inv); } SparseComplexMatrix SparseComplexMatrix::ipermute (const Array<octave_idx_type>& vec) const { return MSparse<Complex>::ipermute (vec); } // other operations bool SparseComplexMatrix::any_element_is_nan (void) const { octave_idx_type nel = nnz (); for (octave_idx_type i = 0; i < nel; i++) { Complex val = data (i); if (octave::math::isnan (val)) return true; } return false; } bool SparseComplexMatrix::any_element_is_inf_or_nan (void) const { octave_idx_type nel = nnz (); for (octave_idx_type i = 0; i < nel; i++) { Complex val = data (i); if (octave::math::isinf (val) || octave::math::isnan (val)) return true; } return false; } // Return true if no elements have imaginary components. bool SparseComplexMatrix::all_elements_are_real (void) const { return mx_inline_all_real (nnz (), data ()); } // Return nonzero if any element of CM has a non-integer real or // imaginary part. Also extract the largest and smallest (real or // imaginary) values and return them in MAX_VAL and MIN_VAL. bool SparseComplexMatrix::all_integers (double& max_val, double& min_val) const { octave_idx_type nel = nnz (); if (nel == 0) return false; max_val = octave::math::real (data (0)); min_val = octave::math::real (data (0)); for (octave_idx_type i = 0; i < nel; i++) { Complex val = data (i); double r_val = val.real (); double i_val = val.imag (); if (r_val > max_val) max_val = r_val; if (i_val > max_val) max_val = i_val; if (r_val < min_val) min_val = r_val; if (i_val < min_val) min_val = i_val; if (octave::math::x_nint (r_val) != r_val || octave::math::x_nint (i_val) != i_val) return false; } return true; } bool SparseComplexMatrix::too_large_for_float (void) const { return test_any (xtoo_large_for_float); } // FIXME: Do these really belong here? Maybe they should be in a base class? SparseBoolMatrix SparseComplexMatrix::all (int dim) const { SPARSE_ALL_OP (dim); } SparseBoolMatrix SparseComplexMatrix::any (int dim) const { SPARSE_ANY_OP (dim); } SparseComplexMatrix SparseComplexMatrix::cumprod (int dim) const { SPARSE_CUMPROD (SparseComplexMatrix, Complex, cumprod); } SparseComplexMatrix SparseComplexMatrix::cumsum (int dim) const { SPARSE_CUMSUM (SparseComplexMatrix, Complex, cumsum); } SparseComplexMatrix SparseComplexMatrix::prod (int dim) const { if ((rows () == 1 && dim == -1) || dim == 1) return transpose ().prod (0).transpose (); else { SPARSE_REDUCTION_OP (SparseComplexMatrix, Complex, *=, (cidx (j+1) - cidx (j) < nr ? 0.0 : 1.0), 1.0); } } SparseComplexMatrix SparseComplexMatrix::sum (int dim) const { SPARSE_REDUCTION_OP (SparseComplexMatrix, Complex, +=, 0.0, 0.0); } SparseComplexMatrix SparseComplexMatrix::sumsq (int dim) const { #define ROW_EXPR \ Complex d = data (i); \ tmp[ridx (i)] += d * conj (d) #define COL_EXPR \ Complex d = data (i); \ tmp[j] += d * conj (d) SPARSE_BASE_REDUCTION_OP (SparseComplexMatrix, Complex, ROW_EXPR, COL_EXPR, 0.0, 0.0); #undef ROW_EXPR #undef COL_EXPR } SparseMatrix SparseComplexMatrix::abs (void) const { octave_idx_type nz = nnz (); octave_idx_type nc = cols (); SparseMatrix retval (rows (), nc, nz); for (octave_idx_type i = 0; i < nc + 1; i++) retval.cidx (i) = cidx (i); for (octave_idx_type i = 0; i < nz; i++) { retval.data (i) = std::abs (data (i)); retval.ridx (i) = ridx (i); } return retval; } SparseComplexMatrix SparseComplexMatrix::diag (octave_idx_type k) const { return MSparse<Complex>::diag (k); } std::ostream& operator << (std::ostream& os, const SparseComplexMatrix& a) { octave_idx_type nc = a.cols (); // add one to the printed indices to go from // zero-based to one-based arrays for (octave_idx_type j = 0; j < nc; j++) { octave_quit (); for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) { os << a.ridx (i) + 1 << " " << j + 1 << " "; octave_write_complex (os, a.data (i)); os << "\n"; } } return os; } std::istream& operator >> (std::istream& is, SparseComplexMatrix& a) { typedef SparseComplexMatrix::element_type elt_type; return read_sparse_matrix<elt_type> (is, a, octave_read_value<Complex>); } SparseComplexMatrix operator * (const SparseComplexMatrix& m, const SparseMatrix& a) { SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, double); } SparseComplexMatrix operator * (const SparseMatrix& m, const SparseComplexMatrix& a) { SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, Complex); } SparseComplexMatrix operator * (const SparseComplexMatrix& m, const SparseComplexMatrix& a) { SPARSE_SPARSE_MUL (SparseComplexMatrix, Complex, Complex); } ComplexMatrix operator * (const ComplexMatrix& m, const SparseMatrix& a) { FULL_SPARSE_MUL (ComplexMatrix, double, Complex (0.,0.)); } ComplexMatrix operator * (const Matrix& m, const SparseComplexMatrix& a) { FULL_SPARSE_MUL (ComplexMatrix, Complex, Complex (0.,0.)); } ComplexMatrix operator * (const ComplexMatrix& m, const SparseComplexMatrix& a) { FULL_SPARSE_MUL (ComplexMatrix, Complex, Complex (0.,0.)); } ComplexMatrix mul_trans (const ComplexMatrix& m, const SparseComplexMatrix& a) { FULL_SPARSE_MUL_TRANS (ComplexMatrix, Complex, Complex (0.,0.), ); } ComplexMatrix mul_herm (const ComplexMatrix& m, const SparseComplexMatrix& a) { FULL_SPARSE_MUL_TRANS (ComplexMatrix, Complex, Complex (0.,0.), conj); } ComplexMatrix operator * (const SparseComplexMatrix& m, const Matrix& a) { SPARSE_FULL_MUL (ComplexMatrix, double, Complex (0.,0.)); } ComplexMatrix operator * (const SparseMatrix& m, const ComplexMatrix& a) { SPARSE_FULL_MUL (ComplexMatrix, Complex, Complex (0.,0.)); } ComplexMatrix operator * (const SparseComplexMatrix& m, const ComplexMatrix& a) { SPARSE_FULL_MUL (ComplexMatrix, Complex, Complex (0.,0.)); } ComplexMatrix trans_mul (const SparseComplexMatrix& m, const ComplexMatrix& a) { SPARSE_FULL_TRANS_MUL (ComplexMatrix, Complex, Complex (0.,0.), ); } ComplexMatrix herm_mul (const SparseComplexMatrix& m, const ComplexMatrix& a) { SPARSE_FULL_TRANS_MUL (ComplexMatrix, Complex, Complex (0.,0.), conj); } // diag * sparse and sparse * diag SparseComplexMatrix operator * (const DiagMatrix& d, const SparseComplexMatrix& a) { return do_mul_dm_sm<SparseComplexMatrix> (d, a); } SparseComplexMatrix operator * (const SparseComplexMatrix& a, const DiagMatrix& d) { return do_mul_sm_dm<SparseComplexMatrix> (a, d); } SparseComplexMatrix operator * (const ComplexDiagMatrix& d, const SparseMatrix& a) { return do_mul_dm_sm<SparseComplexMatrix> (d, a); } SparseComplexMatrix operator * (const SparseMatrix& a, const ComplexDiagMatrix& d) { return do_mul_sm_dm<SparseComplexMatrix> (a, d); } SparseComplexMatrix operator * (const ComplexDiagMatrix& d, const SparseComplexMatrix& a) { return do_mul_dm_sm<SparseComplexMatrix> (d, a); } SparseComplexMatrix operator * (const SparseComplexMatrix& a, const ComplexDiagMatrix& d) { return do_mul_sm_dm<SparseComplexMatrix> (a, d); } SparseComplexMatrix operator + (const ComplexDiagMatrix& d, const SparseMatrix& a) { return do_add_dm_sm<SparseComplexMatrix> (d, a); } SparseComplexMatrix operator + (const DiagMatrix& d, const SparseComplexMatrix& a) { return do_add_dm_sm<SparseComplexMatrix> (d, a); } SparseComplexMatrix operator + (const ComplexDiagMatrix& d, const SparseComplexMatrix& a) { return do_add_dm_sm<SparseComplexMatrix> (d, a); } SparseComplexMatrix operator + (const SparseMatrix& a, const ComplexDiagMatrix& d) { return do_add_sm_dm<SparseComplexMatrix> (a, d); } SparseComplexMatrix operator + (const SparseComplexMatrix& a, const DiagMatrix& d) { return do_add_sm_dm<SparseComplexMatrix> (a, d); } SparseComplexMatrix operator + (const SparseComplexMatrix& a, const ComplexDiagMatrix& d) { return do_add_sm_dm<SparseComplexMatrix> (a, d); } SparseComplexMatrix operator - (const ComplexDiagMatrix& d, const SparseMatrix& a) { return do_sub_dm_sm<SparseComplexMatrix> (d, a); } SparseComplexMatrix operator - (const DiagMatrix& d, const SparseComplexMatrix& a) { return do_sub_dm_sm<SparseComplexMatrix> (d, a); } SparseComplexMatrix operator - (const ComplexDiagMatrix& d, const SparseComplexMatrix& a) { return do_sub_dm_sm<SparseComplexMatrix> (d, a); } SparseComplexMatrix operator - (const SparseMatrix& a, const ComplexDiagMatrix& d) { return do_sub_sm_dm<SparseComplexMatrix> (a, d); } SparseComplexMatrix operator - (const SparseComplexMatrix& a, const DiagMatrix& d) { return do_sub_sm_dm<SparseComplexMatrix> (a, d); } SparseComplexMatrix operator - (const SparseComplexMatrix& a, const ComplexDiagMatrix& d) { return do_sub_sm_dm<SparseComplexMatrix> (a, d); } // perm * sparse and sparse * perm SparseComplexMatrix operator * (const PermMatrix& p, const SparseComplexMatrix& a) { return octinternal_do_mul_pm_sm (p, a); } SparseComplexMatrix operator * (const SparseComplexMatrix& a, const PermMatrix& p) { return octinternal_do_mul_sm_pm (a, p); } // FIXME: it would be nice to share code among the min/max functions below. #define EMPTY_RETURN_CHECK(T) \ if (nr == 0 || nc == 0) \ return T (nr, nc); SparseComplexMatrix min (const Complex& c, const SparseComplexMatrix& m) { SparseComplexMatrix result; octave_idx_type nr = m.rows (); octave_idx_type nc = m.columns (); EMPTY_RETURN_CHECK (SparseComplexMatrix); if (abs (c) == 0.) return SparseComplexMatrix (nr, nc); else { result = SparseComplexMatrix (m); for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) result.data (i) = octave::math::min (c, m.data (i)); } return result; } SparseComplexMatrix min (const SparseComplexMatrix& m, const Complex& c) { return min (c, m); } SparseComplexMatrix min (const SparseComplexMatrix& a, const SparseComplexMatrix& b) { SparseComplexMatrix r; octave_idx_type a_nr = a.rows (); octave_idx_type a_nc = a.cols (); octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); if (a_nr == b_nr && a_nc == b_nc) { r = SparseComplexMatrix (a_nr, a_nc, (a.nnz () + b.nnz ())); octave_idx_type jx = 0; r.cidx (0) = 0; for (octave_idx_type i = 0 ; i < a_nc ; i++) { octave_idx_type ja = a.cidx (i); octave_idx_type ja_max = a.cidx (i+1); bool ja_lt_max = ja < ja_max; octave_idx_type jb = b.cidx (i); octave_idx_type jb_max = b.cidx (i+1); bool jb_lt_max = jb < jb_max; while (ja_lt_max || jb_lt_max) { octave_quit (); if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb)))) { Complex tmp = octave::math::min (a.data (ja), 0.); if (tmp != 0.) { r.ridx (jx) = a.ridx (ja); r.data (jx) = tmp; jx++; } ja++; ja_lt_max= ja < ja_max; } else if ((! ja_lt_max) || (jb_lt_max && (b.ridx (jb) < a.ridx (ja)))) { Complex tmp = octave::math::min (0., b.data (jb)); if (tmp != 0.) { r.ridx (jx) = b.ridx (jb); r.data (jx) = tmp; jx++; } jb++; jb_lt_max= jb < jb_max; } else { Complex tmp = octave::math::min (a.data (ja), b.data (jb)); if (tmp != 0.) { r.data (jx) = tmp; r.ridx (jx) = a.ridx (ja); jx++; } ja++; ja_lt_max= ja < ja_max; jb++; jb_lt_max= jb < jb_max; } } r.cidx (i+1) = jx; } r.maybe_compress (); } else { if (a_nr == 0 || a_nc == 0) r.resize (a_nr, a_nc); else if (b_nr == 0 || b_nc == 0) r.resize (b_nr, b_nc); else octave::err_nonconformant ("min", a_nr, a_nc, b_nr, b_nc); } return r; } SparseComplexMatrix max (const Complex& c, const SparseComplexMatrix& m) { SparseComplexMatrix result; octave_idx_type nr = m.rows (); octave_idx_type nc = m.columns (); EMPTY_RETURN_CHECK (SparseComplexMatrix); // Count the number of nonzero elements if (octave::math::max (c, 0.) != 0.) { result = SparseComplexMatrix (nr, nc, c); for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) result.xdata (m.ridx (i) + j * nr) = octave::math::max (c, m.data (i)); } else result = SparseComplexMatrix (m); return result; } SparseComplexMatrix max (const SparseComplexMatrix& m, const Complex& c) { return max (c, m); } SparseComplexMatrix max (const SparseComplexMatrix& a, const SparseComplexMatrix& b) { SparseComplexMatrix r; octave_idx_type a_nr = a.rows (); octave_idx_type a_nc = a.cols (); octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); if (a_nr == b_nr && a_nc == b_nc) { r = SparseComplexMatrix (a_nr, a_nc, (a.nnz () + b.nnz ())); octave_idx_type jx = 0; r.cidx (0) = 0; for (octave_idx_type i = 0 ; i < a_nc ; i++) { octave_idx_type ja = a.cidx (i); octave_idx_type ja_max = a.cidx (i+1); bool ja_lt_max = ja < ja_max; octave_idx_type jb = b.cidx (i); octave_idx_type jb_max = b.cidx (i+1); bool jb_lt_max = jb < jb_max; while (ja_lt_max || jb_lt_max) { octave_quit (); if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb)))) { Complex tmp = octave::math::max (a.data (ja), 0.); if (tmp != 0.) { r.ridx (jx) = a.ridx (ja); r.data (jx) = tmp; jx++; } ja++; ja_lt_max= ja < ja_max; } else if ((! ja_lt_max) || (jb_lt_max && (b.ridx (jb) < a.ridx (ja)))) { Complex tmp = octave::math::max (0., b.data (jb)); if (tmp != 0.) { r.ridx (jx) = b.ridx (jb); r.data (jx) = tmp; jx++; } jb++; jb_lt_max= jb < jb_max; } else { Complex tmp = octave::math::max (a.data (ja), b.data (jb)); if (tmp != 0.) { r.data (jx) = tmp; r.ridx (jx) = a.ridx (ja); jx++; } ja++; ja_lt_max= ja < ja_max; jb++; jb_lt_max= jb < jb_max; } } r.cidx (i+1) = jx; } r.maybe_compress (); } else { if (a_nr == 0 || a_nc == 0) r.resize (a_nr, a_nc); else if (b_nr == 0 || b_nc == 0) r.resize (b_nr, b_nc); else octave::err_nonconformant ("max", a_nr, a_nc, b_nr, b_nc); } return r; } SPARSE_SMS_CMP_OPS (SparseComplexMatrix, 0.0, real, Complex, 0.0, real) SPARSE_SMS_BOOL_OPS (SparseComplexMatrix, Complex, 0.0) SPARSE_SSM_CMP_OPS (Complex, 0.0, real, SparseComplexMatrix, 0.0, real) SPARSE_SSM_BOOL_OPS (Complex, SparseComplexMatrix, 0.0) SPARSE_SMSM_CMP_OPS (SparseComplexMatrix, 0.0, real, SparseComplexMatrix, 0.0, real) SPARSE_SMSM_BOOL_OPS (SparseComplexMatrix, SparseComplexMatrix, 0.0)