Mercurial > octave
view liboctave/array/Array-base.cc @ 31771:21f9b34eb893
maint: Eliminate "(void)" in C++ function prototypes/declarations.
* mk-opts.pl, external.txi, embedded.cc, make_int.cc, standalone.cc,
standalonebuiltin.cc, BaseControl.cc, BaseControl.h, ButtonControl.cc,
ButtonControl.h, ButtonGroup.cc, ButtonGroup.h, Canvas.cc, Canvas.h,
CheckBoxControl.cc, CheckBoxControl.h, Container.cc, Container.h,
ContextMenu.cc, ContextMenu.h, EditControl.cc, EditControl.h, Figure.cc,
Figure.h, FigureWindow.cc, FigureWindow.h, GLCanvas.cc, GLCanvas.h,
GenericEventNotify.h, KeyMap.cc, ListBoxControl.cc, ListBoxControl.h,
Logger.cc, Logger.h, Menu.cc, Menu.h, MenuContainer.h, Object.cc, Object.h,
ObjectProxy.cc, ObjectProxy.h, Panel.cc, Panel.h, PopupMenuControl.cc,
PopupMenuControl.h, PushButtonControl.cc, PushButtonControl.h, PushTool.cc,
PushTool.h, RadioButtonControl.cc, RadioButtonControl.h, SliderControl.cc,
SliderControl.h, Table.cc, Table.h, TextControl.cc, TextControl.h, TextEdit.h,
ToggleButtonControl.cc, ToggleButtonControl.h, ToggleTool.cc, ToggleTool.h,
ToolBar.cc, ToolBar.h, ToolBarButton.cc, ToolBarButton.h, gl-select.cc,
gl-select.h, qopengl-functions.h, qt-graphics-toolkit.h, qdialog.cpp,
qfontdialog.cpp, qprintdialog_win.cpp, liboctgui-build-info.h,
liboctgui-build-info.in.cc, color-picker.cc, color-picker.h, command-widget.cc,
command-widget.h, community-news.cc, community-news.h, dialog.cc, dialog.h,
documentation-bookmarks.cc, documentation-bookmarks.h,
documentation-dock-widget.cc, documentation-dock-widget.h, documentation.cc,
documentation.h, dw-main-window.cc, dw-main-window.h,
external-editor-interface.cc, external-editor-interface.h,
files-dock-widget.cc, files-dock-widget.h, find-files-dialog.cc,
find-files-dialog.h, find-files-model.cc, find-files-model.h,
gui-preferences.cc, gui-preferences.h, gui-settings.cc, gui-settings.h,
history-dock-widget.cc, history-dock-widget.h, interpreter-qobject.cc,
interpreter-qobject.h, file-editor-interface.h, file-editor-tab.cc,
file-editor-tab.h, file-editor.cc, file-editor.h, find-dialog.cc,
find-dialog.h, marker.cc, marker.h, octave-qscintilla.cc, octave-qscintilla.h,
octave-txt-lexer.cc, octave-txt-lexer.h, main-window.cc, main-window.h,
news-reader.cc, news-reader.h, octave-dock-widget.cc, octave-dock-widget.h,
octave-qobject.cc, octave-qobject.h, qt-application.cc, qt-application.h,
qt-interpreter-events.cc, qt-interpreter-events.h, release-notes.cc,
release-notes.h, set-path-dialog.cc, set-path-dialog.h, set-path-model.cc,
set-path-model.h, settings-dialog.cc, settings-dialog.h,
shortcuts-tree-widget.cc, shortcuts-tree-widget.h, tab-bar.cc, tab-bar.h,
terminal-dock-widget.cc, terminal-dock-widget.h, variable-editor-model.cc,
variable-editor-model.h, variable-editor.cc, variable-editor.h,
welcome-wizard.cc, welcome-wizard.h, workspace-model.cc, workspace-model.h,
workspace-view.cc, workspace-view.h, build-env.h, Cell.cc, Cell.h,
__contourc__.cc, __magick_read__.cc, auto-shlib.cc, auto-shlib.h,
base-text-renderer.h, bsxfun.cc, c-file-ptr-stream.cc, c-file-ptr-stream.h,
call-stack.cc, call-stack.h, debug.cc, defaults.cc, defaults.h, defun.cc,
display.cc, display.h, dynamic-ld.cc, dynamic-ld.h, environment.cc,
environment.h, error.cc, error.h, errwarn.cc, errwarn.h, event-manager.cc,
event-manager.h, event-queue.cc, event-queue.h, fcn-info.cc, fcn-info.h,
ft-text-renderer.cc, ft-text-renderer.h, genprops.awk, gh-manager.cc,
gh-manager.h, gl-render.cc, gl-render.h, gl2ps-print.cc, graphics-toolkit.h,
graphics.cc, graphics.in.h, gtk-manager.cc, gtk-manager.h, help.cc, help.h,
hook-fcn.h, input.cc, input.h, interpreter-private.cc, interpreter-private.h,
interpreter.cc, interpreter.h, jsondecode.cc, latex-text-renderer.cc,
latex-text-renderer.h, load-path.cc, load-path.h, load-save.cc, load-save.h,
ls-hdf5.cc, ls-hdf5.h, mxarray.h, oct-errno.h, oct-errno.in.cc, oct-fstrm.cc,
oct-fstrm.h, oct-handle.h, oct-hist.cc, oct-hist.h, oct-iostrm.cc,
oct-iostrm.h, oct-map.cc, oct-map.h, oct-opengl.h, oct-prcstrm.cc,
oct-prcstrm.h, oct-procbuf.cc, oct-procbuf.h, oct-process.h, oct-stdstrm.h,
oct-stream.cc, oct-stream.h, oct-strstrm.cc, oct-strstrm.h,
oct-tex-lexer.in.ll, pager.cc, pager.h, pr-flt-fmt.cc, pr-flt-fmt.h,
pr-output.cc, pr-output.h, procstream.cc, procstream.h, settings.cc,
settings.h, sighandlers.cc, sighandlers.h, stack-frame.cc, stack-frame.h,
svd.cc, syminfo.cc, syminfo.h, symrec.cc, symrec.h, symscope.cc, symscope.h,
symtab.cc, symtab.h, sysdep.cc, sysdep.h, text-engine.cc, text-engine.h,
text-renderer.cc, text-renderer.h, toplev.cc, url-handle-manager.cc,
url-handle-manager.h, variables.cc, xpow.cc, __init_fltk__.cc,
__init_gnuplot__.cc, __ode15__.cc, audiodevinfo.cc, gzip.cc,
liboctinterp-build-info.h, liboctinterp-build-info.in.cc,
mk-build-env-features.sh, mk-builtins.pl, cdef-class.cc, cdef-class.h,
cdef-manager.h, cdef-method.cc, cdef-method.h, cdef-object.cc, cdef-object.h,
cdef-package.cc, cdef-package.h, cdef-property.cc, cdef-property.h,
cdef-utils.cc, cdef-utils.h, ov-base-diag.cc, ov-base-diag.h, ov-base-int.cc,
ov-base-int.h, ov-base-mat.cc, ov-base-mat.h, ov-base-scalar.cc,
ov-base-scalar.h, ov-base-sparse.cc, ov-base-sparse.h, ov-base.cc, ov-base.h,
ov-bool-mat.cc, ov-bool-mat.h, ov-bool-sparse.cc, ov-bool-sparse.h, ov-bool.cc,
ov-bool.h, ov-builtin.cc, ov-builtin.h, ov-cell.cc, ov-cell.h, ov-ch-mat.cc,
ov-ch-mat.h, ov-class.cc, ov-class.h, ov-classdef.cc, ov-classdef.h,
ov-colon.h, ov-complex.cc, ov-complex.h, ov-cs-list.h, ov-cx-diag.cc,
ov-cx-diag.h, ov-cx-mat.cc, ov-cx-mat.h, ov-cx-sparse.cc, ov-cx-sparse.h,
ov-dld-fcn.cc, ov-dld-fcn.h, ov-fcn-handle.cc, ov-fcn-handle.h, ov-fcn.cc,
ov-fcn.h, ov-float.cc, ov-float.h, ov-flt-complex.cc, ov-flt-complex.h,
ov-flt-cx-diag.cc, ov-flt-cx-diag.h, ov-flt-cx-mat.cc, ov-flt-cx-mat.h,
ov-flt-re-diag.cc, ov-flt-re-diag.h, ov-flt-re-mat.cc, ov-flt-re-mat.h,
ov-intx.h, ov-java.cc, ov-java.h, ov-lazy-idx.cc, ov-lazy-idx.h,
ov-legacy-range.cc, ov-legacy-range.h, ov-magic-int.cc, ov-magic-int.h,
ov-mex-fcn.cc, ov-mex-fcn.h, ov-null-mat.cc, ov-null-mat.h, ov-oncleanup.cc,
ov-oncleanup.h, ov-perm.cc, ov-perm.h, ov-range.cc, ov-range.h, ov-re-diag.cc,
ov-re-diag.h, ov-re-mat.cc, ov-re-mat.h, ov-re-sparse.cc, ov-re-sparse.h,
ov-scalar.cc, ov-scalar.h, ov-str-mat.cc, ov-str-mat.h, ov-struct.cc,
ov-struct.h, ov-typeinfo.cc, ov-typeinfo.h, ov-usr-fcn.cc, ov-usr-fcn.h, ov.cc,
ov.h, ovl.cc, ovl.h, octave.cc, octave.h, anon-fcn-validator.h, bp-table.cc,
bp-table.h, comment-list.cc, comment-list.h, filepos.h, lex.h, lex.ll,
oct-lvalue.cc, oct-lvalue.h, oct-parse.yy, parse.h, profiler.cc, profiler.h,
pt-anon-scopes.h, pt-arg-list.cc, pt-arg-list.h, pt-args-block.cc,
pt-args-block.h, pt-array-list.cc, pt-array-list.h, pt-assign.cc, pt-assign.h,
pt-binop.cc, pt-binop.h, pt-bp.h, pt-cbinop.h, pt-cell.h, pt-check.h,
pt-classdef.cc, pt-classdef.h, pt-cmd.h, pt-colon.h, pt-const.h, pt-decl.cc,
pt-decl.h, pt-eval.cc, pt-eval.h, pt-except.cc, pt-except.h, pt-exp.cc,
pt-exp.h, pt-fcn-handle.cc, pt-fcn-handle.h, pt-id.cc, pt-id.h, pt-idx.cc,
pt-idx.h, pt-jump.h, pt-loop.cc, pt-loop.h, pt-mat.h, pt-misc.cc, pt-misc.h,
pt-pr-code.cc, pt-pr-code.h, pt-select.cc, pt-select.h, pt-spmd.cc, pt-spmd.h,
pt-stmt.cc, pt-stmt.h, pt-tm-const.cc, pt-tm-const.h, pt-unop.cc, pt-unop.h,
pt-walk.h, pt.cc, pt.h, token.cc, token.h, usage.h, Array-base.cc, Array.h,
CColVector.cc, CColVector.h, CDiagMatrix.cc, CDiagMatrix.h, CMatrix.cc,
CMatrix.h, CNDArray.cc, CNDArray.h, CRowVector.cc, CRowVector.h, CSparse.cc,
CSparse.h, DiagArray2.cc, DiagArray2.h, MArray.cc, MArray.h, MDiagArray2.h,
MSparse.h, MatrixType.cc, MatrixType.h, PermMatrix.cc, PermMatrix.h, Range.cc,
Range.h, Sparse-b.cc, Sparse.cc, Sparse.h, boolMatrix.cc, boolMatrix.h,
boolNDArray.cc, boolNDArray.h, boolSparse.cc, boolSparse.h, chMatrix.h,
chNDArray.h, dColVector.cc, dColVector.h, dDiagMatrix.cc, dDiagMatrix.h,
dMatrix.cc, dMatrix.h, dNDArray.cc, dNDArray.h, dRowVector.cc, dRowVector.h,
dSparse.cc, dSparse.h, dim-vector.cc, dim-vector.h, fCColVector.cc,
fCColVector.h, fCDiagMatrix.cc, fCDiagMatrix.h, fCMatrix.cc, fCMatrix.h,
fCNDArray.cc, fCNDArray.h, fCRowVector.cc, fCRowVector.h, fColVector.cc,
fColVector.h, fDiagMatrix.cc, fDiagMatrix.h, fMatrix.cc, fMatrix.h,
fNDArray.cc, fNDArray.h, fRowVector.cc, fRowVector.h, idx-vector.cc,
idx-vector.h, intNDArray.cc, intNDArray.h, liboctave-build-info.h,
liboctave-build-info.in.cc, CollocWt.cc, CollocWt.h, DAE.h, DAEFunc.h, DAERT.h,
DAERTFunc.h, DASPK.cc, DASPK.h, DASRT.cc, DASRT.h, DASSL.cc, DASSL.h, DET.h,
EIG.h, LSODE.cc, LSODE.h, ODE.h, ODEFunc.h, ODES.h, ODESFunc.h, Quad.h,
aepbalance.cc, aepbalance.h, base-dae.h, base-de.h, chol.cc, chol.h,
eigs-base.cc, fEIG.h, gepbalance.h, gsvd.cc, gsvd.h, hess.h, lu.cc, lu.h,
oct-fftw.cc, oct-fftw.h, oct-rand.cc, oct-rand.h, oct-spparms.cc,
oct-spparms.h, qr.cc, qr.h, qrp.cc, qrp.h, randmtzig.cc, randmtzig.h, schur.h,
sparse-chol.cc, sparse-chol.h, sparse-lu.cc, sparse-lu.h, sparse-qr.cc,
sparse-qr.h, svd.cc, svd.h, child-list.cc, child-list.h, dir-ops.cc, dir-ops.h,
file-ops.cc, file-ops.h, file-stat.cc, file-stat.h, lo-sysdep.cc, lo-sysdep.h,
lo-sysinfo.cc, lo-sysinfo.h, mach-info.cc, mach-info.h, oct-env.cc, oct-env.h,
oct-group.cc, oct-group.h, oct-password.cc, oct-password.h, oct-syscalls.cc,
oct-syscalls.h, oct-time.cc, oct-time.h, oct-uname.cc, oct-uname.h,
action-container.h, base-list.h, caseless-str.h, cmd-edit.cc, cmd-edit.h,
cmd-hist.cc, cmd-hist.h, data-conv.cc, file-info.h, glob-match.cc,
glob-match.h, kpse.cc, kpse.h, lo-array-errwarn.cc, lo-array-errwarn.h,
lo-hash.cc, lo-hash.h, lo-ieee.cc, lo-ieee.h, lo-regexp.cc, lo-regexp.h,
oct-inttypes.cc, oct-inttypes.h, oct-mutex.cc, oct-mutex.h, oct-refcount.h,
oct-shlib.cc, oct-shlib.h, oct-sort.cc, oct-sort.h, oct-sparse.cc,
octave-preserve-stream-state.h, pathsearch.cc, pathsearch.h, quit.cc,
singleton-cleanup.cc, singleton-cleanup.h, str-vec.cc, str-vec.h,
unwind-prot.cc, unwind-prot.h, url-transfer.cc, url-transfer.h, version.cc,
version.in.h, cxx-signal-helpers.cc, acinclude.m4, main-cli.cc, main-gui.cc,
main.in.cc, mkoctfile.in.cc, octave-build-info.h, octave-build-info.in.cc,
octave-config.in.cc, octave-svgconvert.cc, shared-fcns.h:
maint: Eliminate "(void)" in C++ function prototypes/declarations.
author | Rik <rik@octave.org> |
---|---|
date | Tue, 24 Jan 2023 17:19:44 -0800 |
parents | 597f3ee61a48 |
children | febd82d1a8de |
line wrap: on
line source
//////////////////////////////////////////////////////////////////////// // // Copyright (C) 1993-2023 The Octave Project Developers // // See the file COPYRIGHT.md in the top-level directory of this // distribution or <https://octave.org/copyright/>. // // 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 // <https://www.gnu.org/licenses/>. // //////////////////////////////////////////////////////////////////////// // Do not include this file directly to instantiate the Array<T> template // class in code outside liboctave. Include "Array-oct.cc" or "Array.cc" // instead. // This file should not include config.h. It is only included in other // C++ source files that should have included config.h before including // this file. #include <ostream> #include "Array-util.h" #include "Array.h" #include "lo-mappers.h" #include "oct-locbuf.h" // One dimensional array class. Handles the reference counting for // all the derived classes. template <typename T, typename Alloc> typename Array<T, Alloc>::ArrayRep * Array<T, Alloc>::nil_rep () { static ArrayRep nr; return &nr; } // This is similar to the template for containers but specialized for Array. // Note we can't specialize a member without also specializing the class. template <typename T, typename Alloc> Array<T, Alloc>::Array (const Array<T, Alloc>& a, const dim_vector& dv) : m_dimensions (dv), m_rep (a.m_rep), m_slice_data (a.m_slice_data), m_slice_len (a.m_slice_len) { bool invalid_size = false; octave_idx_type new_numel = 0; try { // Safe numel may throw a bad_alloc error because the new // dimensions are too large to allocate. Trap that here so // we can issue a more helpful diagnostic that is consistent // with other size mismatch failures. new_numel = m_dimensions.safe_numel (); } catch (const std::bad_alloc&) { invalid_size = true; } if (invalid_size || new_numel != a.numel ()) { std::string dimensions_str = a.m_dimensions.str (); std::string new_dims_str = m_dimensions.str (); (*current_liboctave_error_handler) ("reshape: can't reshape %s array to %s array", dimensions_str.c_str (), new_dims_str.c_str ()); } // This goes here because if an exception is thrown by the above, // destructor will be never called. m_rep->m_count++; m_dimensions.chop_trailing_singletons (); } template <typename T, typename Alloc> void Array<T, Alloc>::fill (const T& val) { if (m_rep->m_count > 1) { --m_rep->m_count; m_rep = new ArrayRep (numel (), val); m_slice_data = m_rep->m_data; } else std::fill_n (m_slice_data, m_slice_len, val); } template <typename T, typename Alloc> void Array<T, Alloc>::clear () { if (--m_rep->m_count == 0) delete m_rep; m_rep = nil_rep (); m_rep->m_count++; m_slice_data = m_rep->m_data; m_slice_len = m_rep->m_len; m_dimensions = dim_vector (); } template <typename T, typename Alloc> void Array<T, Alloc>::clear (const dim_vector& dv) { if (--m_rep->m_count == 0) delete m_rep; m_rep = new ArrayRep (dv.safe_numel ()); m_slice_data = m_rep->m_data; m_slice_len = m_rep->m_len; m_dimensions = dv; m_dimensions.chop_trailing_singletons (); } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::squeeze () const { Array<T, Alloc> retval = *this; if (ndims () > 2) { bool dims_changed = false; dim_vector new_dimensions = m_dimensions; int k = 0; for (int i = 0; i < ndims (); i++) { if (m_dimensions(i) == 1) dims_changed = true; else new_dimensions(k++) = m_dimensions(i); } if (dims_changed) { switch (k) { case 0: new_dimensions = dim_vector (1, 1); break; case 1: { octave_idx_type tmp = new_dimensions(0); new_dimensions.resize (2); new_dimensions(0) = tmp; new_dimensions(1) = 1; } break; default: new_dimensions.resize (k); break; } } retval = Array<T, Alloc> (*this, new_dimensions); } return retval; } template <typename T, typename Alloc> octave_idx_type Array<T, Alloc>::compute_index (octave_idx_type i, octave_idx_type j) const { return ::compute_index (i, j, m_dimensions); } template <typename T, typename Alloc> octave_idx_type Array<T, Alloc>::compute_index (octave_idx_type i, octave_idx_type j, octave_idx_type k) const { return ::compute_index (i, j, k, m_dimensions); } template <typename T, typename Alloc> octave_idx_type Array<T, Alloc>::compute_index (const Array<octave_idx_type>& ra_idx) const { return ::compute_index (ra_idx, m_dimensions); } template <typename T, typename Alloc> T& Array<T, Alloc>::checkelem (octave_idx_type n) { // Do checks directly to avoid recomputing m_slice_len. if (n < 0) octave::err_invalid_index (n); if (n >= m_slice_len) octave::err_index_out_of_range (1, 1, n+1, m_slice_len, m_dimensions); return elem (n); } template <typename T, typename Alloc> T& Array<T, Alloc>::checkelem (octave_idx_type i, octave_idx_type j) { return elem (compute_index (i, j)); } template <typename T, typename Alloc> T& Array<T, Alloc>::checkelem (octave_idx_type i, octave_idx_type j, octave_idx_type k) { return elem (compute_index (i, j, k)); } template <typename T, typename Alloc> T& Array<T, Alloc>::checkelem (const Array<octave_idx_type>& ra_idx) { return elem (compute_index (ra_idx)); } template <typename T, typename Alloc> typename Array<T, Alloc>::crefT Array<T, Alloc>::checkelem (octave_idx_type n) const { // Do checks directly to avoid recomputing m_slice_len. if (n < 0) octave::err_invalid_index (n); if (n >= m_slice_len) octave::err_index_out_of_range (1, 1, n+1, m_slice_len, m_dimensions); return elem (n); } template <typename T, typename Alloc> typename Array<T, Alloc>::crefT Array<T, Alloc>::checkelem (octave_idx_type i, octave_idx_type j) const { return elem (compute_index (i, j)); } template <typename T, typename Alloc> typename Array<T, Alloc>::crefT Array<T, Alloc>::checkelem (octave_idx_type i, octave_idx_type j, octave_idx_type k) const { return elem (compute_index (i, j, k)); } template <typename T, typename Alloc> typename Array<T, Alloc>::crefT Array<T, Alloc>::checkelem (const Array<octave_idx_type>& ra_idx) const { return elem (compute_index (ra_idx)); } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::column (octave_idx_type k) const { octave_idx_type r = m_dimensions(0); return Array<T, Alloc> (*this, dim_vector (r, 1), k*r, k*r + r); } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::page (octave_idx_type k) const { octave_idx_type r = m_dimensions(0); octave_idx_type c = m_dimensions(1); octave_idx_type p = r*c; return Array<T, Alloc> (*this, dim_vector (r, c), k*p, k*p + p); } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::linear_slice (octave_idx_type lo, octave_idx_type up) const { if (up < lo) up = lo; return Array<T, Alloc> (*this, dim_vector (up - lo, 1), lo, up); } // Helper class for multi-d dimension permuting (generalized transpose). class rec_permute_helper { public: rec_permute_helper (const dim_vector& dv, const Array<octave_idx_type>& perm) : m_n (dv.ndims ()), m_top (0), m_dim (new octave_idx_type [2*m_n]), m_stride (m_dim + m_n), m_use_blk (false) { assert (m_n == perm.numel ()); // Get cumulative dimensions. OCTAVE_LOCAL_BUFFER (octave_idx_type, cdim, m_n+1); cdim[0] = 1; for (int i = 1; i < m_n+1; i++) cdim[i] = cdim[i-1] * dv(i-1); // Setup the permuted strides. for (int k = 0; k < m_n; k++) { int kk = perm(k); m_dim[k] = dv(kk); m_stride[k] = cdim[kk]; } // Reduce contiguous runs. for (int k = 1; k < m_n; k++) { if (m_stride[k] == m_stride[m_top]*m_dim[m_top]) m_dim[m_top] *= m_dim[k]; else { m_top++; m_dim[m_top] = m_dim[k]; m_stride[m_top] = m_stride[k]; } } // Determine whether we can use block transposes. m_use_blk = m_top >= 1 && m_stride[1] == 1 && m_stride[0] == m_dim[1]; } // No copying! rec_permute_helper (const rec_permute_helper&) = delete; rec_permute_helper& operator = (const rec_permute_helper&) = delete; ~rec_permute_helper () { delete [] m_dim; } template <typename T> void permute (const T *src, T *dest) const { do_permute (src, dest, m_top); } // Helper method for fast blocked transpose. template <typename T> static T * blk_trans (const T *src, T *dest, octave_idx_type nr, octave_idx_type nc) { static const octave_idx_type m = 8; OCTAVE_LOCAL_BUFFER (T, blk, m*m); for (octave_idx_type kr = 0; kr < nr; kr += m) for (octave_idx_type kc = 0; kc < nc; kc += m) { octave_idx_type lr = std::min (m, nr - kr); octave_idx_type lc = std::min (m, nc - kc); if (lr == m && lc == m) { const T *ss = src + kc * nr + kr; for (octave_idx_type j = 0; j < m; j++) for (octave_idx_type i = 0; i < m; i++) blk[j*m+i] = ss[j*nr + i]; T *dd = dest + kr * nc + kc; for (octave_idx_type j = 0; j < m; j++) for (octave_idx_type i = 0; i < m; i++) dd[j*nc+i] = blk[i*m+j]; } else { const T *ss = src + kc * nr + kr; for (octave_idx_type j = 0; j < lc; j++) for (octave_idx_type i = 0; i < lr; i++) blk[j*m+i] = ss[j*nr + i]; T *dd = dest + kr * nc + kc; for (octave_idx_type j = 0; j < lr; j++) for (octave_idx_type i = 0; i < lc; i++) dd[j*nc+i] = blk[i*m+j]; } } return dest + nr*nc; } private: // Recursive N-D generalized transpose template <typename T> T * do_permute (const T *src, T *dest, int lev) const { if (lev == 0) { octave_idx_type step = m_stride[0]; octave_idx_type len = m_dim[0]; if (step == 1) { std::copy_n (src, len, dest); dest += len; } else { for (octave_idx_type i = 0, j = 0; i < len; i++, j += step) dest[i] = src[j]; dest += len; } } else if (m_use_blk && lev == 1) dest = blk_trans (src, dest, m_dim[1], m_dim[0]); else { octave_idx_type step = m_stride[lev]; octave_idx_type len = m_dim[lev]; for (octave_idx_type i = 0, j = 0; i < len; i++, j+= step) dest = do_permute (src + i * step, dest, lev-1); } return dest; } //-------- // STRIDE occupies the last half of the space allocated for dim to // avoid a double allocation. int m_n; int m_top; octave_idx_type *m_dim; octave_idx_type *m_stride; bool m_use_blk; }; template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::permute (const Array<octave_idx_type>& perm_vec_arg, bool inv) const { Array<T, Alloc> retval; Array<octave_idx_type> perm_vec = perm_vec_arg; dim_vector dv = dims (); int perm_vec_len = perm_vec_arg.numel (); if (perm_vec_len < dv.ndims ()) (*current_liboctave_error_handler) ("%s: invalid permutation vector", inv ? "ipermute" : "permute"); dim_vector dv_new = dim_vector::alloc (perm_vec_len); // Append singleton dimensions as needed. dv.resize (perm_vec_len, 1); // Need this array to check for identical elements in permutation array. OCTAVE_LOCAL_BUFFER_INIT (bool, checked, perm_vec_len, false); bool identity = true; // Find dimension vector of permuted array. for (int i = 0; i < perm_vec_len; i++) { octave_idx_type perm_elt = perm_vec.elem (i); if (perm_elt >= perm_vec_len || perm_elt < 0) (*current_liboctave_error_handler) ("%s: permutation vector contains an invalid element", inv ? "ipermute" : "permute"); if (checked[perm_elt]) (*current_liboctave_error_handler) ("%s: permutation vector cannot contain identical elements", inv ? "ipermute" : "permute"); else { checked[perm_elt] = true; identity = identity && perm_elt == i; } } if (identity) return *this; if (inv) { for (int i = 0; i < perm_vec_len; i++) perm_vec(perm_vec_arg(i)) = i; } for (int i = 0; i < perm_vec_len; i++) dv_new(i) = dv(perm_vec(i)); retval = Array<T, Alloc> (dv_new); if (numel () > 0) { rec_permute_helper rh (dv, perm_vec); rh.permute (data (), retval.fortran_vec ()); } return retval; } // Helper class for multi-d index reduction and recursive // indexing/indexed assignment. Rationale: we could avoid recursion // using a state machine instead. However, using recursion is much // more amenable to possible parallelization in the future. // Also, the recursion solution is cleaner and more understandable. class rec_index_helper { public: rec_index_helper (const dim_vector& dv, const Array<octave::idx_vector>& ia) : m_n (ia.numel ()), m_top (0), m_dim (new octave_idx_type [2*m_n]), m_cdim (m_dim + m_n), m_idx (new octave::idx_vector [m_n]) { assert (m_n > 0 && (dv.ndims () == std::max (m_n, 2))); m_dim[0] = dv(0); m_cdim[0] = 1; m_idx[0] = ia(0); for (int i = 1; i < m_n; i++) { // Try reduction... if (m_idx[m_top].maybe_reduce (m_dim[m_top], ia(i), dv(i))) { // Reduction successful, fold dimensions. m_dim[m_top] *= dv(i); } else { // Unsuccessful, store index & cumulative m_dim. m_top++; m_idx[m_top] = ia(i); m_dim[m_top] = dv(i); m_cdim[m_top] = m_cdim[m_top-1] * m_dim[m_top-1]; } } } // No copying! rec_index_helper (const rec_index_helper&) = delete; rec_index_helper& operator = (const rec_index_helper&) = delete; ~rec_index_helper () { delete [] m_idx; delete [] m_dim; } template <typename T> void index (const T *src, T *dest) const { do_index (src, dest, m_top); } template <typename T> void assign (const T *src, T *dest) const { do_assign (src, dest, m_top); } template <typename T> void fill (const T& val, T *dest) const { do_fill (val, dest, m_top); } bool is_cont_range (octave_idx_type& l, octave_idx_type& u) const { return m_top == 0 && m_idx[0].is_cont_range (m_dim[0], l, u); } private: // Recursive N-D indexing template <typename T> T * do_index (const T *src, T *dest, int lev) const { if (lev == 0) dest += m_idx[0].index (src, m_dim[0], dest); else { octave_idx_type nn = m_idx[lev].length (m_dim[lev]); octave_idx_type d = m_cdim[lev]; for (octave_idx_type i = 0; i < nn; i++) dest = do_index (src + d*m_idx[lev].xelem (i), dest, lev-1); } return dest; } // Recursive N-D indexed assignment template <typename T> const T * do_assign (const T *src, T *dest, int lev) const { if (lev == 0) src += m_idx[0].assign (src, m_dim[0], dest); else { octave_idx_type nn = m_idx[lev].length (m_dim[lev]); octave_idx_type d = m_cdim[lev]; for (octave_idx_type i = 0; i < nn; i++) src = do_assign (src, dest + d*m_idx[lev].xelem (i), lev-1); } return src; } // Recursive N-D indexed assignment template <typename T> void do_fill (const T& val, T *dest, int lev) const { if (lev == 0) m_idx[0].fill (val, m_dim[0], dest); else { octave_idx_type nn = m_idx[lev].length (m_dim[lev]); octave_idx_type d = m_cdim[lev]; for (octave_idx_type i = 0; i < nn; i++) do_fill (val, dest + d*m_idx[lev].xelem (i), lev-1); } } //-------- // CDIM occupies the last half of the space allocated for dim to // avoid a double allocation. int m_n; int m_top; octave_idx_type *m_dim; octave_idx_type *m_cdim; octave::idx_vector *m_idx; }; // Helper class for multi-d recursive resizing // This handles resize () in an efficient manner, touching memory only // once (apart from reinitialization) class rec_resize_helper { public: rec_resize_helper (const dim_vector& ndv, const dim_vector& odv) : m_cext (nullptr), m_sext (nullptr), m_dext (nullptr), m_n (0) { int l = ndv.ndims (); assert (odv.ndims () == l); octave_idx_type ld = 1; int i = 0; for (; i < l-1 && ndv(i) == odv(i); i++) ld *= ndv(i); m_n = l - i; m_cext = new octave_idx_type [3*m_n]; // Trick to avoid three allocations m_sext = m_cext + m_n; m_dext = m_sext + m_n; octave_idx_type sld = ld; octave_idx_type dld = ld; for (int j = 0; j < m_n; j++) { m_cext[j] = std::min (ndv(i+j), odv(i+j)); m_sext[j] = sld *= odv(i+j); m_dext[j] = dld *= ndv(i+j); } m_cext[0] *= ld; } // No copying! rec_resize_helper (const rec_resize_helper&) = delete; rec_resize_helper& operator = (const rec_resize_helper&) = delete; ~rec_resize_helper () { delete [] m_cext; } template <typename T> void resize_fill (const T *src, T *dest, const T& rfv) const { do_resize_fill (src, dest, rfv, m_n-1); } private: // recursive resizing template <typename T> void do_resize_fill (const T *src, T *dest, const T& rfv, int lev) const { if (lev == 0) { std::copy_n (src, m_cext[0], dest); std::fill_n (dest + m_cext[0], m_dext[0] - m_cext[0], rfv); } else { octave_idx_type sd, dd, k; sd = m_sext[lev-1]; dd = m_dext[lev-1]; for (k = 0; k < m_cext[lev]; k++) do_resize_fill (src + k * sd, dest + k * dd, rfv, lev - 1); std::fill_n (dest + k * dd, m_dext[lev] - k * dd, rfv); } } //-------- octave_idx_type *m_cext; octave_idx_type *m_sext; octave_idx_type *m_dext; int m_n; }; template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::index (const octave::idx_vector& i) const { // Colon: // // object | index | result orientation // ---------+----------+------------------- // anything | colon | column vector // // Numeric array or logical mask: // // Note that logical mask indices are always transformed to vectors // before we see them. // // object | index | result orientation // ---------+----------+------------------- // vector | vector | indexed object // | other | same size as index // ---------+----------+------------------- // array | anything | same size as index octave_idx_type n = numel (); Array<T, Alloc> retval; if (i.is_colon ()) { // A(:) produces a shallow copy as a column vector. retval = Array<T, Alloc> (*this, dim_vector (n, 1)); } else { if (i.extent (n) != n) octave::err_index_out_of_range (1, 1, i.extent (n), n, m_dimensions); dim_vector result_dims = i.orig_dimensions (); octave_idx_type idx_len = i.length (); if (n != 1 && is_nd_vector () && idx_len != 1 && result_dims.is_nd_vector ()) { // Indexed object and index are both vectors. Set result size // and orientation as above. dim_vector dv = dims (); result_dims = dv.make_nd_vector (idx_len); } octave_idx_type l, u; if (idx_len != 0 && i.is_cont_range (n, l, u)) // If suitable, produce a shallow slice. retval = Array<T, Alloc> (*this, result_dims, l, u); else { // Don't use resize here to avoid useless initialization for POD // types. retval = Array<T, Alloc> (result_dims); if (idx_len != 0) i.index (data (), n, retval.fortran_vec ()); } } return retval; } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::index (const octave::idx_vector& i, const octave::idx_vector& j) const { // Get dimensions, allowing Fortran indexing in the 2nd dim. dim_vector dv = m_dimensions.redim (2); octave_idx_type r = dv(0); octave_idx_type c = dv(1); Array<T, Alloc> retval; if (i.is_colon () && j.is_colon ()) { // A(:,:) produces a shallow copy. retval = Array<T, Alloc> (*this, dv); } else { if (i.extent (r) != r) octave::err_index_out_of_range (2, 1, i.extent (r), r, m_dimensions); // throws if (j.extent (c) != c) octave::err_index_out_of_range (2, 2, j.extent (c), c, m_dimensions); // throws octave_idx_type n = numel (); octave_idx_type il = i.length (r); octave_idx_type jl = j.length (c); octave::idx_vector ii (i); if (ii.maybe_reduce (r, j, c)) { octave_idx_type l, u; if (ii.length () > 0 && ii.is_cont_range (n, l, u)) // If suitable, produce a shallow slice. retval = Array<T, Alloc> (*this, dim_vector (il, jl), l, u); else { // Don't use resize to avoid useless initialization for POD types. retval = Array<T, Alloc> (dim_vector (il, jl)); ii.index (data (), n, retval.fortran_vec ()); } } else { // Don't use resize to avoid useless initialization for POD types. retval = Array<T, Alloc> (dim_vector (il, jl)); const T *src = data (); T *dest = retval.fortran_vec (); for (octave_idx_type k = 0; k < jl; k++) dest += i.index (src + r * j.xelem (k), r, dest); } } return retval; } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::index (const Array<octave::idx_vector>& ia) const { int ial = ia.numel (); Array<T, Alloc> retval; // FIXME: is this dispatching necessary? if (ial == 1) retval = index (ia(0)); else if (ial == 2) retval = index (ia(0), ia(1)); else if (ial > 0) { // Get dimensions, allowing Fortran indexing in the last dim. dim_vector dv = m_dimensions.redim (ial); // Check for out of bounds conditions. bool all_colons = true; for (int i = 0; i < ial; i++) { if (ia(i).extent (dv(i)) != dv(i)) octave::err_index_out_of_range (ial, i+1, ia(i).extent (dv(i)), dv(i), m_dimensions); // throws all_colons = all_colons && ia(i).is_colon (); } if (all_colons) { // A(:,:,...,:) produces a shallow copy. dv.chop_trailing_singletons (); retval = Array<T, Alloc> (*this, dv); } else { // Form result dimensions. dim_vector rdv = dim_vector::alloc (ial); for (int i = 0; i < ial; i++) rdv(i) = ia(i).length (dv(i)); rdv.chop_trailing_singletons (); // Prepare for recursive indexing rec_index_helper rh (dv, ia); octave_idx_type l, u; if (rh.is_cont_range (l, u)) // If suitable, produce a shallow slice. retval = Array<T, Alloc> (*this, rdv, l, u); else { // Don't use resize to avoid useless initialization for POD types. retval = Array<T, Alloc> (rdv); // Do it. rh.index (data (), retval.fortran_vec ()); } } } return retval; } // The default fill value. Override if you want a different one. template <typename T, typename Alloc> T Array<T, Alloc>::resize_fill_value () const { static T zero = T (); return zero; } // Yes, we could do resize using index & assign. However, that would // possibly involve a lot more memory traffic than we actually need. template <typename T, typename Alloc> void Array<T, Alloc>::resize1 (octave_idx_type n, const T& rfv) { if (n < 0 || ndims () != 2) octave::err_invalid_resize (); dim_vector dv; // This is driven by Matlab's behavior of giving a *row* vector // on some out-of-bounds assignments. Specifically, Matlab // allows a(i) with out-of-bouds i when a is either of 0x0, 1x0, // 1x1, 0xN, and gives a row vector in all cases (yes, even the // last one, search me why). Giving a column vector would make // much more sense (given the way trailing singleton dims are // treated). bool invalid = false; if (rows () == 0 || rows () == 1) dv = dim_vector (1, n); else if (columns () == 1) dv = dim_vector (n, 1); else invalid = true; if (invalid) octave::err_invalid_resize (); octave_idx_type nx = numel (); if (n == nx - 1 && n > 0) { // Stack "pop" operation. if (m_rep->m_count == 1) m_slice_data[m_slice_len-1] = T (); m_slice_len--; m_dimensions = dv; } else if (n == nx + 1 && nx > 0) { // Stack "push" operation. if (m_rep->m_count == 1 && m_slice_data + m_slice_len < m_rep->m_data + m_rep->m_len) { m_slice_data[m_slice_len++] = rfv; m_dimensions = dv; } else { static const octave_idx_type max_stack_chunk = 1024; octave_idx_type nn = n + std::min (nx, max_stack_chunk); Array<T, Alloc> tmp (Array<T, Alloc> (dim_vector (nn, 1)), dv, 0, n); T *dest = tmp.fortran_vec (); std::copy_n (data (), nx, dest); dest[nx] = rfv; *this = tmp; } } else if (n != nx) { Array<T, Alloc> tmp = Array<T, Alloc> (dv); T *dest = tmp.fortran_vec (); octave_idx_type n0 = std::min (n, nx); octave_idx_type n1 = n - n0; std::copy_n (data (), n0, dest); std::fill_n (dest + n0, n1, rfv); *this = tmp; } } template <typename T, typename Alloc> void Array<T, Alloc>::resize2 (octave_idx_type r, octave_idx_type c, const T& rfv) { if (r < 0 || c < 0 || ndims () != 2) octave::err_invalid_resize (); octave_idx_type rx = rows (); octave_idx_type cx = columns (); if (r != rx || c != cx) { Array<T, Alloc> tmp = Array<T, Alloc> (dim_vector (r, c)); T *dest = tmp.fortran_vec (); octave_idx_type r0 = std::min (r, rx); octave_idx_type r1 = r - r0; octave_idx_type c0 = std::min (c, cx); octave_idx_type c1 = c - c0; const T *src = data (); if (r == rx) { std::copy_n (src, r * c0, dest); dest += r * c0; } else { for (octave_idx_type k = 0; k < c0; k++) { std::copy_n (src, r0, dest); src += rx; dest += r0; std::fill_n (dest, r1, rfv); dest += r1; } } std::fill_n (dest, r * c1, rfv); *this = tmp; } } template <typename T, typename Alloc> void Array<T, Alloc>::resize (const dim_vector& dv, const T& rfv) { int dvl = dv.ndims (); if (dvl == 2) resize2 (dv(0), dv(1), rfv); else if (m_dimensions != dv) { if (m_dimensions.ndims () > dvl || dv.any_neg ()) octave::err_invalid_resize (); Array<T, Alloc> tmp (dv); // Prepare for recursive resizing. rec_resize_helper rh (dv, m_dimensions.redim (dvl)); // Do it. rh.resize_fill (data (), tmp.fortran_vec (), rfv); *this = tmp; } } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::index (const octave::idx_vector& i, bool resize_ok, const T& rfv) const { Array<T, Alloc> tmp = *this; if (resize_ok) { octave_idx_type n = numel (); octave_idx_type nx = i.extent (n); if (n != nx) { if (i.is_scalar ()) return Array<T, Alloc> (dim_vector (1, 1), rfv); else tmp.resize1 (nx, rfv); } if (tmp.numel () != nx) return Array<T, Alloc> (); } return tmp.index (i); } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::index (const octave::idx_vector& i, const octave::idx_vector& j, bool resize_ok, const T& rfv) const { Array<T, Alloc> tmp = *this; if (resize_ok) { dim_vector dv = m_dimensions.redim (2); octave_idx_type r = dv(0); octave_idx_type c = dv(1); octave_idx_type rx = i.extent (r); octave_idx_type cx = j.extent (c); if (r != rx || c != cx) { if (i.is_scalar () && j.is_scalar ()) return Array<T, Alloc> (dim_vector (1, 1), rfv); else tmp.resize2 (rx, cx, rfv); } if (tmp.rows () != rx || tmp.columns () != cx) return Array<T, Alloc> (); } return tmp.index (i, j); } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::index (const Array<octave::idx_vector>& ia, bool resize_ok, const T& rfv) const { Array<T, Alloc> tmp = *this; if (resize_ok) { int ial = ia.numel (); dim_vector dv = m_dimensions.redim (ial); dim_vector dvx = dim_vector::alloc (ial); for (int i = 0; i < ial; i++) dvx(i) = ia(i).extent (dv(i)); if (! (dvx == dv)) { bool all_scalars = true; for (int i = 0; i < ial; i++) all_scalars = all_scalars && ia(i).is_scalar (); if (all_scalars) return Array<T, Alloc> (dim_vector (1, 1), rfv); else tmp.resize (dvx, rfv); if (tmp.m_dimensions != dvx) return Array<T, Alloc> (); } } return tmp.index (ia); } template <typename T, typename Alloc> void Array<T, Alloc>::assign (const octave::idx_vector& i, const Array<T, Alloc>& rhs, const T& rfv) { octave_idx_type n = numel (); octave_idx_type rhl = rhs.numel (); if (rhl != 1 && i.length (n) != rhl) octave::err_nonconformant ("=", dim_vector(i.length(n), 1), rhs.dims()); octave_idx_type nx = i.extent (n); bool colon = i.is_colon_equiv (nx); // Try to resize first if necessary. if (nx != n) { // Optimize case A = []; A(1:n) = X with A empty. if (m_dimensions.zero_by_zero () && colon) { if (rhl == 1) *this = Array<T, Alloc> (dim_vector (1, nx), rhs(0)); else *this = Array<T, Alloc> (rhs, dim_vector (1, nx)); return; } resize1 (nx, rfv); n = numel (); } if (colon) { // A(:) = X makes a full fill or a shallow copy. if (rhl == 1) fill (rhs(0)); else *this = rhs.reshape (m_dimensions); } else { if (rhl == 1) i.fill (rhs(0), n, fortran_vec ()); else i.assign (rhs.data (), n, fortran_vec ()); } } // Assignment to a 2-dimensional array template <typename T, typename Alloc> void Array<T, Alloc>::assign (const octave::idx_vector& i, const octave::idx_vector& j, const Array<T, Alloc>& rhs, const T& rfv) { bool initial_dims_all_zero = m_dimensions.all_zero (); // Get RHS extents, discarding singletons. dim_vector rhdv = rhs.dims (); // Get LHS extents, allowing Fortran indexing in the second dim. dim_vector dv = m_dimensions.redim (2); // Check for out-of-bounds and form resizing m_dimensions. dim_vector rdv; // In the special when all dimensions are zero, colons are allowed // to inquire the shape of RHS. The rules are more obscure, so we // solve that elsewhere. if (initial_dims_all_zero) rdv = zero_dims_inquire (i, j, rhdv); else { rdv(0) = i.extent (dv(0)); rdv(1) = j.extent (dv(1)); } bool isfill = rhs.numel () == 1; octave_idx_type il = i.length (rdv(0)); octave_idx_type jl = j.length (rdv(1)); rhdv.chop_all_singletons (); bool match = (isfill || (rhdv.ndims () == 2 && il == rhdv(0) && jl == rhdv(1))); match = match || (il == 1 && jl == rhdv(0) && rhdv(1) == 1); if (match) { bool all_colons = (i.is_colon_equiv (rdv(0)) && j.is_colon_equiv (rdv(1))); // Resize if requested. if (rdv != dv) { // Optimize case A = []; A(1:m, 1:n) = X if (dv.zero_by_zero () && all_colons) { if (isfill) *this = Array<T, Alloc> (rdv, rhs(0)); else *this = Array<T, Alloc> (rhs, rdv); return; } resize (rdv, rfv); dv = m_dimensions; } if (all_colons) { // A(:,:) = X makes a full fill or a shallow copy if (isfill) fill (rhs(0)); else *this = rhs.reshape (m_dimensions); } else { // The actual work. octave_idx_type n = numel (); octave_idx_type r = dv(0); octave_idx_type c = dv(1); octave::idx_vector ii (i); const T *src = rhs.data (); T *dest = fortran_vec (); // Try reduction first. if (ii.maybe_reduce (r, j, c)) { if (isfill) ii.fill (*src, n, dest); else ii.assign (src, n, dest); } else { if (isfill) { for (octave_idx_type k = 0; k < jl; k++) i.fill (*src, r, dest + r * j.xelem (k)); } else { for (octave_idx_type k = 0; k < jl; k++) src += i.assign (src, r, dest + r * j.xelem (k)); } } } } // any empty RHS can be assigned to an empty LHS else if ((il != 0 && jl != 0) || (rhdv(0) != 0 && rhdv(1) != 0)) octave::err_nonconformant ("=", il, jl, rhs.dim1 (), rhs.dim2 ()); } // Assignment to a multi-dimensional array template <typename T, typename Alloc> void Array<T, Alloc>::assign (const Array<octave::idx_vector>& ia, const Array<T, Alloc>& rhs, const T& rfv) { int ial = ia.numel (); // FIXME: is this dispatching necessary / desirable? if (ial == 1) assign (ia(0), rhs, rfv); else if (ial == 2) assign (ia(0), ia(1), rhs, rfv); else if (ial > 0) { bool initial_dims_all_zero = m_dimensions.all_zero (); // Get RHS extents, discarding singletons. dim_vector rhdv = rhs.dims (); // Get LHS extents, allowing Fortran indexing in the second dim. dim_vector dv = m_dimensions.redim (ial); // Get the extents forced by indexing. dim_vector rdv; // In the special when all dimensions are zero, colons are // allowed to inquire the shape of RHS. The rules are more // obscure, so we solve that elsewhere. if (initial_dims_all_zero) rdv = zero_dims_inquire (ia, rhdv); else { rdv = dim_vector::alloc (ial); for (int i = 0; i < ial; i++) rdv(i) = ia(i).extent (dv(i)); } // Check whether LHS and RHS match, up to singleton dims. bool match = true; bool all_colons = true; bool isfill = rhs.numel () == 1; rhdv.chop_all_singletons (); int j = 0; int rhdvl = rhdv.ndims (); for (int i = 0; i < ial; i++) { all_colons = all_colons && ia(i).is_colon_equiv (rdv(i)); octave_idx_type l = ia(i).length (rdv(i)); if (l == 1) continue; match = match && j < rhdvl && l == rhdv(j++); } match = match && (j == rhdvl || rhdv(j) == 1); match = match || isfill; if (match) { // Resize first if necessary. if (rdv != dv) { // Optimize case A = []; A(1:m, 1:n) = X if (dv.zero_by_zero () && all_colons) { rdv.chop_trailing_singletons (); if (isfill) *this = Array<T, Alloc> (rdv, rhs(0)); else *this = Array<T, Alloc> (rhs, rdv); return; } resize (rdv, rfv); dv = rdv; } if (all_colons) { // A(:,:,...,:) = X makes a full fill or a shallow copy. if (isfill) fill (rhs(0)); else *this = rhs.reshape (m_dimensions); } else { // Do the actual work. // Prepare for recursive indexing rec_index_helper rh (dv, ia); // Do it. if (isfill) rh.fill (rhs(0), fortran_vec ()); else rh.assign (rhs.data (), fortran_vec ()); } } else { // dimension mismatch, unless LHS and RHS both empty bool lhsempty, rhsempty; lhsempty = rhsempty = false; dim_vector lhs_dv = dim_vector::alloc (ial); for (int i = 0; i < ial; i++) { octave_idx_type l = ia(i).length (rdv(i)); lhs_dv(i) = l; lhsempty = lhsempty || (l == 0); rhsempty = rhsempty || (rhdv(j++) == 0); } if (! lhsempty || ! rhsempty) { lhs_dv.chop_trailing_singletons (); octave::err_nonconformant ("=", lhs_dv, rhdv); } } } } /* %!shared a %! a = [1 2; 3 4]; %!error <op1 is 1x2, op2 is 1x3> a(1,:) = [1 2 3] %!error <op1 is 2x1, op2 is 1x3> a(:,1) = [1 2 3] %!error <op1 is 2x2, op2 is 2x1> a(1:2,2:3) = [1;2] */ template <typename T, typename Alloc> void Array<T, Alloc>::delete_elements (const octave::idx_vector& i) { octave_idx_type n = numel (); if (i.is_colon ()) { *this = Array<T, Alloc> (); } else if (i.length (n) != 0) { if (i.extent (n) != n) octave::err_del_index_out_of_range (true, i.extent (n), n); octave_idx_type l, u; bool col_vec = ndims () == 2 && columns () == 1 && rows () != 1; if (i.is_scalar () && i(0) == n-1 && m_dimensions.isvector ()) { // Stack "pop" operation. resize1 (n-1); } else if (i.is_cont_range (n, l, u)) { // Special case deleting a contiguous range. octave_idx_type m = n + l - u; Array<T, Alloc> tmp (dim_vector (col_vec ? m : 1, ! col_vec ? m : 1)); const T *src = data (); T *dest = tmp.fortran_vec (); std::copy_n (src, l, dest); std::copy (src + u, src + n, dest + l); *this = tmp; } else { // Use index. *this = index (i.complement (n)); } } } template <typename T, typename Alloc> void Array<T, Alloc>::delete_elements (int dim, const octave::idx_vector& i) { if (dim < 0 || dim >= ndims ()) (*current_liboctave_error_handler) ("invalid dimension in delete_elements"); octave_idx_type n = m_dimensions(dim); if (i.is_colon ()) { *this = Array<T, Alloc> (); } else if (i.length (n) != 0) { if (i.extent (n) != n) octave::err_del_index_out_of_range (false, i.extent (n), n); octave_idx_type l, u; if (i.is_cont_range (n, l, u)) { // Special case deleting a contiguous range. octave_idx_type nd = n + l - u; octave_idx_type dl = 1; octave_idx_type du = 1; dim_vector rdv = m_dimensions; rdv(dim) = nd; for (int k = 0; k < dim; k++) dl *= m_dimensions(k); for (int k = dim + 1; k < ndims (); k++) du *= m_dimensions(k); // Special case deleting a contiguous range. Array<T, Alloc> tmp = Array<T, Alloc> (rdv); const T *src = data (); T *dest = tmp.fortran_vec (); l *= dl; u *= dl; n *= dl; for (octave_idx_type k = 0; k < du; k++) { std::copy_n (src, l, dest); dest += l; std::copy (src + u, src + n, dest); dest += n - u; src += n; } *this = tmp; } else { // Use index. Array<octave::idx_vector> ia (dim_vector (ndims (), 1), octave::idx_vector::colon); ia (dim) = i.complement (n); *this = index (ia); } } } template <typename T, typename Alloc> void Array<T, Alloc>::delete_elements (const Array<octave::idx_vector>& ia) { int ial = ia.numel (); if (ial == 1) delete_elements (ia(0)); else { int k, dim = -1; for (k = 0; k < ial; k++) { if (! ia(k).is_colon ()) { if (dim < 0) dim = k; else break; } } if (dim < 0) { dim_vector dv = m_dimensions; dv(0) = 0; *this = Array<T, Alloc> (dv); } else if (k == ial) { delete_elements (dim, ia(dim)); } else { // Allow the null assignment to succeed if it won't change // anything because the indices reference an empty slice, // provided that there is at most one non-colon (or // equivalent) index. So, we still have the requirement of // deleting a slice, but it is OK if the slice is empty. // For compatibility with Matlab, stop checking once we see // more than one non-colon index or an empty index. Matlab // considers "[]" to be an empty index but not "false". We // accept both. bool empty_assignment = false; int num_non_colon_indices = 0; int nd = ndims (); for (int i = 0; i < ial; i++) { octave_idx_type dim_len = (i >= nd ? 1 : m_dimensions(i)); if (ia(i).length (dim_len) == 0) { empty_assignment = true; break; } if (! ia(i).is_colon_equiv (dim_len)) { num_non_colon_indices++; if (num_non_colon_indices == 2) break; } } if (! empty_assignment) (*current_liboctave_error_handler) ("a null assignment can only have one non-colon index"); } } } template <typename T, typename Alloc> Array<T, Alloc>& Array<T, Alloc>::insert (const Array<T, Alloc>& a, octave_idx_type r, octave_idx_type c) { octave::idx_vector i (r, r + a.rows ()); octave::idx_vector j (c, c + a.columns ()); if (ndims () == 2 && a.ndims () == 2) assign (i, j, a); else { Array<octave::idx_vector> idx (dim_vector (a.ndims (), 1)); idx(0) = i; idx(1) = j; for (int k = 2; k < a.ndims (); k++) idx(k) = octave::idx_vector (0, a.m_dimensions(k)); assign (idx, a); } return *this; } template <typename T, typename Alloc> Array<T, Alloc>& Array<T, Alloc>::insert (const Array<T, Alloc>& a, const Array<octave_idx_type>& ra_idx) { octave_idx_type n = ra_idx.numel (); Array<octave::idx_vector> idx (dim_vector (n, 1)); const dim_vector dva = a.dims ().redim (n); for (octave_idx_type k = 0; k < n; k++) idx(k) = octave::idx_vector (ra_idx(k), ra_idx(k) + dva(k)); assign (idx, a); return *this; } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::transpose () const { assert (ndims () == 2); octave_idx_type nr = dim1 (); octave_idx_type nc = dim2 (); if (nr >= 8 && nc >= 8) { Array<T, Alloc> result (dim_vector (nc, nr)); // Reuse the implementation used for permuting. rec_permute_helper::blk_trans (data (), result.fortran_vec (), nr, nc); return result; } else if (nr > 1 && nc > 1) { Array<T, Alloc> result (dim_vector (nc, nr)); for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = 0; i < nr; i++) result.xelem (j, i) = xelem (i, j); return result; } else { // Fast transpose for vectors and empty matrices. return Array<T, Alloc> (*this, dim_vector (nc, nr)); } } template <typename T> static T no_op_fcn (const T& x) { return x; } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::hermitian (T (*fcn) (const T&)) const { assert (ndims () == 2); if (! fcn) fcn = no_op_fcn<T>; octave_idx_type nr = dim1 (); octave_idx_type nc = dim2 (); if (nr >= 8 && nc >= 8) { Array<T, Alloc> result (dim_vector (nc, nr)); // Blocked transpose to attempt to avoid cache misses. T buf[64]; octave_idx_type jj; for (jj = 0; jj < (nc - 8 + 1); jj += 8) { octave_idx_type ii; for (ii = 0; ii < (nr - 8 + 1); ii += 8) { // Copy to buffer for (octave_idx_type j = jj, k = 0, idxj = jj * nr; j < jj + 8; j++, idxj += nr) for (octave_idx_type i = ii; i < ii + 8; i++) buf[k++] = xelem (i + idxj); // Copy from buffer for (octave_idx_type i = ii, idxi = ii * nc; i < ii + 8; i++, idxi += nc) for (octave_idx_type j = jj, k = i - ii; j < jj + 8; j++, k+=8) result.xelem (j + idxi) = fcn (buf[k]); } if (ii < nr) for (octave_idx_type j = jj; j < jj + 8; j++) for (octave_idx_type i = ii; i < nr; i++) result.xelem (j, i) = fcn (xelem (i, j)); } for (octave_idx_type j = jj; j < nc; j++) for (octave_idx_type i = 0; i < nr; i++) result.xelem (j, i) = fcn (xelem (i, j)); return result; } else { Array<T, Alloc> result (dim_vector (nc, nr)); for (octave_idx_type j = 0; j < nc; j++) for (octave_idx_type i = 0; i < nr; i++) result.xelem (j, i) = fcn (xelem (i, j)); return result; } } /* ## Transpose tests for matrices of the tile size and plus or minus a row ## and with four tiles. %!shared m7, mt7, m8, mt8, m9, mt9 %! m7 = reshape (1 : 7*8, 8, 7); %! mt7 = [1:8; 9:16; 17:24; 25:32; 33:40; 41:48; 49:56]; %! m8 = reshape (1 : 8*8, 8, 8); %! mt8 = mt8 = [mt7; 57:64]; %! m9 = reshape (1 : 9*8, 8, 9); %! mt9 = [mt8; 65:72]; %!assert (m7', mt7) %!assert ((1i*m7).', 1i * mt7) %!assert ((1i*m7)', conj (1i * mt7)) %!assert (m8', mt8) %!assert ((1i*m8).', 1i * mt8) %!assert ((1i*m8)', conj (1i * mt8)) %!assert (m9', mt9) %!assert ((1i*m9).', 1i * mt9) %!assert ((1i*m9)', conj (1i * mt9)) %!assert ([m7, m8; m7, m8]', [mt7, mt7; mt8, mt8]) %!assert ((1i*[m7, m8; m7, m8]).', 1i * [mt7, mt7; mt8, mt8]) %!assert ((1i*[m7, m8; m7, m8])', conj (1i * [mt7, mt7; mt8, mt8])) %!assert ([m8, m8; m8, m8]', [mt8, mt8; mt8, mt8]) %!assert ((1i*[m8, m8; m8, m8]).', 1i * [mt8, mt8; mt8, mt8]) %!assert ((1i*[m8, m8; m8, m8])', conj (1i * [mt8, mt8; mt8, mt8])) %!assert ([m9, m8; m9, m8]', [mt9, mt9; mt8, mt8]) %!assert ((1i*[m9, m8; m9, m8]).', 1i * [mt9, mt9; mt8, mt8]) %!assert ((1i*[m9, m8; m9, m8])', conj (1i * [mt9, mt9; mt8, mt8])) */ template <typename T, typename Alloc> T * Array<T, Alloc>::fortran_vec () { make_unique (); return m_slice_data; } // Non-real types don't have NaNs. template <typename T> inline bool sort_isnan (typename ref_param<T>::type) { return false; } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::sort (int dim, sortmode mode) const { if (dim < 0) (*current_liboctave_error_handler) ("sort: invalid dimension"); Array<T, Alloc> m (dims ()); dim_vector dv = m.dims (); if (m.numel () < 1) return m; if (dim >= dv.ndims ()) dv.resize (dim+1, 1); octave_idx_type ns = dv(dim); octave_idx_type iter = dv.numel () / ns; octave_idx_type stride = 1; for (int i = 0; i < dim; i++) stride *= dv(i); T *v = m.fortran_vec (); const T *ov = data (); octave_sort<T> lsort; if (mode != UNSORTED) lsort.set_compare (mode); else return m; if (stride == 1) { // Special case along first dimension avoids gather/scatter AND directly // sorts into destination buffer for an 11% performance boost. for (octave_idx_type j = 0; j < iter; j++) { // Copy and partition out NaNs. // No need to special case integer types <T> from floating point // types <T> to avoid sort_isnan() test as it makes no discernible // performance impact. octave_idx_type kl = 0; octave_idx_type ku = ns; for (octave_idx_type i = 0; i < ns; i++) { T tmp = ov[i]; if (sort_isnan<T> (tmp)) v[--ku] = tmp; else v[kl++] = tmp; } // sort. lsort.sort (v, kl); if (ku < ns) { // NaNs are in reverse order std::reverse (v + ku, v + ns); if (mode == DESCENDING) std::rotate (v, v + ku, v + ns); } v += ns; ov += ns; } } else { OCTAVE_LOCAL_BUFFER (T, buf, ns); for (octave_idx_type j = 0; j < iter; j++) { octave_idx_type offset = j; octave_idx_type n_strides = j / stride; offset += n_strides * stride * (ns - 1); // gather and partition out NaNs. octave_idx_type kl = 0; octave_idx_type ku = ns; for (octave_idx_type i = 0; i < ns; i++) { T tmp = ov[i*stride + offset]; if (sort_isnan<T> (tmp)) buf[--ku] = tmp; else buf[kl++] = tmp; } // sort. lsort.sort (buf, kl); if (ku < ns) { // NaNs are in reverse order std::reverse (buf + ku, buf + ns); if (mode == DESCENDING) std::rotate (buf, buf + ku, buf + ns); } // scatter. for (octave_idx_type i = 0; i < ns; i++) v[i*stride + offset] = buf[i]; } } return m; } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::sort (Array<octave_idx_type>& sidx, int dim, sortmode mode) const { if (dim < 0 || dim >= ndims ()) (*current_liboctave_error_handler) ("sort: invalid dimension"); Array<T, Alloc> m (dims ()); dim_vector dv = m.dims (); if (m.numel () < 1) { sidx = Array<octave_idx_type> (dv); return m; } octave_idx_type ns = dv(dim); octave_idx_type iter = dv.numel () / ns; octave_idx_type stride = 1; for (int i = 0; i < dim; i++) stride *= dv(i); T *v = m.fortran_vec (); const T *ov = data (); octave_sort<T> lsort; sidx = Array<octave_idx_type> (dv); octave_idx_type *vi = sidx.fortran_vec (); if (mode != UNSORTED) lsort.set_compare (mode); else return m; if (stride == 1) { // Special case for dim 1 avoids gather/scatter for performance boost. // See comments in Array::sort (dim, mode). for (octave_idx_type j = 0; j < iter; j++) { // copy and partition out NaNs. octave_idx_type kl = 0; octave_idx_type ku = ns; for (octave_idx_type i = 0; i < ns; i++) { T tmp = ov[i]; if (sort_isnan<T> (tmp)) { --ku; v[ku] = tmp; vi[ku] = i; } else { v[kl] = tmp; vi[kl] = i; kl++; } } // sort. lsort.sort (v, vi, kl); if (ku < ns) { // NaNs are in reverse order std::reverse (v + ku, v + ns); std::reverse (vi + ku, vi + ns); if (mode == DESCENDING) { std::rotate (v, v + ku, v + ns); std::rotate (vi, vi + ku, vi + ns); } } v += ns; vi += ns; ov += ns; } } else { OCTAVE_LOCAL_BUFFER (T, buf, ns); OCTAVE_LOCAL_BUFFER (octave_idx_type, bufi, ns); for (octave_idx_type j = 0; j < iter; j++) { octave_idx_type offset = j; octave_idx_type n_strides = j / stride; offset += n_strides * stride * (ns - 1); // gather and partition out NaNs. octave_idx_type kl = 0; octave_idx_type ku = ns; for (octave_idx_type i = 0; i < ns; i++) { T tmp = ov[i*stride + offset]; if (sort_isnan<T> (tmp)) { --ku; buf[ku] = tmp; bufi[ku] = i; } else { buf[kl] = tmp; bufi[kl] = i; kl++; } } // sort. lsort.sort (buf, bufi, kl); if (ku < ns) { // NaNs are in reverse order std::reverse (buf + ku, buf + ns); std::reverse (bufi + ku, bufi + ns); if (mode == DESCENDING) { std::rotate (buf, buf + ku, buf + ns); std::rotate (bufi, bufi + ku, bufi + ns); } } // scatter. for (octave_idx_type i = 0; i < ns; i++) v[i*stride + offset] = buf[i]; for (octave_idx_type i = 0; i < ns; i++) vi[i*stride + offset] = bufi[i]; } } return m; } template <typename T, typename Alloc> typename Array<T, Alloc>::compare_fcn_type safe_comparator (sortmode mode, const Array<T, Alloc>& /* a */, bool /* allow_chk */) { if (mode == ASCENDING) return octave_sort<T>::ascending_compare; else if (mode == DESCENDING) return octave_sort<T>::descending_compare; else return nullptr; } template <typename T, typename Alloc> sortmode Array<T, Alloc>::issorted (sortmode mode) const { octave_sort<T> lsort; octave_idx_type n = numel (); if (n <= 1) return (mode == UNSORTED) ? ASCENDING : mode; if (mode == UNSORTED) { // Auto-detect mode. compare_fcn_type compare = safe_comparator (ASCENDING, *this, false); if (compare (elem (n-1), elem (0))) mode = DESCENDING; else mode = ASCENDING; } lsort.set_compare (safe_comparator (mode, *this, false)); if (! lsort.issorted (data (), n)) mode = UNSORTED; return mode; } template <typename T, typename Alloc> Array<octave_idx_type> Array<T, Alloc>::sort_rows_idx (sortmode mode) const { Array<octave_idx_type> idx; octave_sort<T> lsort (safe_comparator (mode, *this, true)); octave_idx_type r = rows (); octave_idx_type c = cols (); idx = Array<octave_idx_type> (dim_vector (r, 1)); lsort.sort_rows (data (), idx.fortran_vec (), r, c); return idx; } template <typename T, typename Alloc> sortmode Array<T, Alloc>::is_sorted_rows (sortmode mode) const { octave_sort<T> lsort; octave_idx_type r = rows (); octave_idx_type c = cols (); if (r <= 1 || c == 0) return (mode == UNSORTED) ? ASCENDING : mode; if (mode == UNSORTED) { // Auto-detect mode. compare_fcn_type compare = safe_comparator (ASCENDING, *this, false); octave_idx_type i; for (i = 0; i < cols (); i++) { T l = elem (0, i); T u = elem (rows () - 1, i); if (compare (l, u)) { if (mode == DESCENDING) { mode = UNSORTED; break; } else mode = ASCENDING; } else if (compare (u, l)) { if (mode == ASCENDING) { mode = UNSORTED; break; } else mode = DESCENDING; } } if (mode == UNSORTED && i == cols ()) mode = ASCENDING; } if (mode != UNSORTED) { lsort.set_compare (safe_comparator (mode, *this, false)); if (! lsort.is_sorted_rows (data (), r, c)) mode = UNSORTED; } return mode; } // Do a binary lookup in a sorted array. template <typename T, typename Alloc> octave_idx_type Array<T, Alloc>::lookup (const T& value, sortmode mode) const { octave_idx_type n = numel (); octave_sort<T> lsort; if (mode == UNSORTED) { // auto-detect mode if (n > 1 && lsort.descending_compare (elem (0), elem (n-1))) mode = DESCENDING; else mode = ASCENDING; } lsort.set_compare (mode); return lsort.lookup (data (), n, value); } template <typename T, typename Alloc> Array<octave_idx_type> Array<T, Alloc>::lookup (const Array<T, Alloc>& values, sortmode mode) const { octave_idx_type n = numel (); octave_idx_type nval = values.numel (); octave_sort<T> lsort; Array<octave_idx_type> idx (values.dims ()); if (mode == UNSORTED) { // auto-detect mode if (n > 1 && lsort.descending_compare (elem (0), elem (n-1))) mode = DESCENDING; else mode = ASCENDING; } lsort.set_compare (mode); // This determines the split ratio between the O(M*log2(N)) and O(M+N) // algorithms. static const double ratio = 1.0; sortmode vmode = UNSORTED; // Attempt the O(M+N) algorithm if M is large enough. if (nval > ratio * n / octave::math::log2 (n + 1.0)) { vmode = values.issorted (); // The table must not contain a NaN. if ((vmode == ASCENDING && sort_isnan<T> (values(nval-1))) || (vmode == DESCENDING && sort_isnan<T> (values(0)))) vmode = UNSORTED; } if (vmode != UNSORTED) lsort.lookup_sorted (data (), n, values.data (), nval, idx.fortran_vec (), vmode != mode); else lsort.lookup (data (), n, values.data (), nval, idx.fortran_vec ()); return idx; } template <typename T, typename Alloc> octave_idx_type Array<T, Alloc>::nnz () const { const T *src = data (); octave_idx_type nel = numel (); octave_idx_type retval = 0; const T zero = T (); for (octave_idx_type i = 0; i < nel; i++) if (src[i] != zero) retval++; return retval; } template <typename T, typename Alloc> Array<octave_idx_type> Array<T, Alloc>::find (octave_idx_type n, bool backward) const { Array<octave_idx_type> retval; const T *src = data (); octave_idx_type nel = numel (); const T zero = T (); if (n < 0 || n >= nel) { // We want all elements, which means we'll almost surely need // to resize. So count first, then allocate array of exact size. octave_idx_type cnt = 0; for (octave_idx_type i = 0; i < nel; i++) cnt += src[i] != zero; retval.clear (cnt, 1); octave_idx_type *dest = retval.fortran_vec (); for (octave_idx_type i = 0; i < nel; i++) if (src[i] != zero) *dest++ = i; } else { // We want a fixed max number of elements, usually small. So be // optimistic, alloc the array in advance, and then resize if // needed. retval.clear (n, 1); if (backward) { // Do the search as a series of successive single-element searches. octave_idx_type k = 0; octave_idx_type l = nel - 1; for (; k < n; k++) { for (; l >= 0 && src[l] == zero; l--) ; if (l >= 0) retval(k) = l--; else break; } if (k < n) retval.resize2 (k, 1); octave_idx_type *rdata = retval.fortran_vec (); std::reverse (rdata, rdata + k); } else { // Do the search as a series of successive single-element searches. octave_idx_type k = 0; octave_idx_type l = 0; for (; k < n; k++) { for (; l != nel && src[l] == zero; l++) ; if (l != nel) retval(k) = l++; else break; } if (k < n) retval.resize2 (k, 1); } } // Fixup return dimensions, for Matlab compatibility. // find (zeros (0,0)) -> zeros (0,0) // find (zeros (1,0)) -> zeros (1,0) // find (zeros (0,1)) -> zeros (0,1) // find (zeros (0,X)) -> zeros (0,1) // find (zeros (1,1)) -> zeros (0,0) !!!! WHY? // find (zeros (0,1,0)) -> zeros (0,0) // find (zeros (0,1,0,1)) -> zeros (0,0) etc if ((numel () == 1 && retval.isempty ()) || (rows () == 0 && dims ().numel (1) == 0)) retval.m_dimensions = dim_vector (); else if (rows () == 1 && ndims () == 2) retval.m_dimensions = dim_vector (1, retval.numel ()); return retval; } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::nth_element (const octave::idx_vector& n, int dim) const { if (dim < 0) (*current_liboctave_error_handler) ("nth_element: invalid dimension"); dim_vector dv = dims (); if (dim >= dv.ndims ()) dv.resize (dim+1, 1); octave_idx_type ns = dv(dim); octave_idx_type nn = n.length (ns); dv(dim) = std::min (nn, ns); dv.chop_trailing_singletons (); dim = std::min (dv.ndims (), static_cast<octave_idx_type> (dim)); Array<T, Alloc> m (dv); if (m.isempty ()) return m; sortmode mode = UNSORTED; octave_idx_type lo = 0; switch (n.idx_class ()) { case octave::idx_vector::class_scalar: mode = ASCENDING; lo = n(0); break; case octave::idx_vector::class_range: { octave_idx_type inc = n.increment (); if (inc == 1) { mode = ASCENDING; lo = n(0); } else if (inc == -1) { mode = DESCENDING; lo = ns - 1 - n(0); } } break; case octave::idx_vector::class_vector: // This case resolves bug #51329, a fallback to allow the given index // to be a sequential vector instead of the typical scalar or range if (n(1) - n(0) == 1) { mode = ASCENDING; lo = n(0); } else if (n(1) - n(0) == -1) { mode = DESCENDING; lo = ns - 1 - n(0); } // Ensure that the vector is actually an arithmetic contiguous sequence for (octave_idx_type i = 2; i < n.length () && mode != UNSORTED; i++) if ((mode == ASCENDING && n(i) - n(i-1) != 1) || (mode == DESCENDING && n(i) - n(i-1) != -1)) mode = UNSORTED; break; default: break; } if (mode == UNSORTED) (*current_liboctave_error_handler) ("nth_element: n must be a scalar or a contiguous range"); octave_idx_type up = lo + nn; if (lo < 0 || up > ns) (*current_liboctave_error_handler) ("nth_element: invalid element index"); octave_idx_type iter = numel () / ns; octave_idx_type stride = 1; for (int i = 0; i < dim; i++) stride *= dv(i); T *v = m.fortran_vec (); const T *ov = data (); OCTAVE_LOCAL_BUFFER (T, buf, ns); octave_sort<T> lsort; lsort.set_compare (mode); for (octave_idx_type j = 0; j < iter; j++) { octave_idx_type kl = 0; octave_idx_type ku = ns; if (stride == 1) { // copy without NaNs. for (octave_idx_type i = 0; i < ns; i++) { T tmp = ov[i]; if (sort_isnan<T> (tmp)) buf[--ku] = tmp; else buf[kl++] = tmp; } ov += ns; } else { octave_idx_type offset = j % stride; // copy without NaNs. for (octave_idx_type i = 0; i < ns; i++) { T tmp = ov[offset + i*stride]; if (sort_isnan<T> (tmp)) buf[--ku] = tmp; else buf[kl++] = tmp; } if (offset == stride-1) ov += ns*stride; } if (ku == ns) lsort.nth_element (buf, ns, lo, up); else if (mode == ASCENDING) lsort.nth_element (buf, ku, lo, std::min (ku, up)); else { octave_idx_type nnan = ns - ku; octave_idx_type zero = 0; lsort.nth_element (buf, ku, std::max (lo - nnan, zero), std::max (up - nnan, zero)); std::rotate (buf, buf + ku, buf + ns); } if (stride == 1) { for (octave_idx_type i = 0; i < nn; i++) v[i] = buf[lo + i]; v += nn; } else { octave_idx_type offset = j % stride; for (octave_idx_type i = 0; i < nn; i++) v[offset + stride * i] = buf[lo + i]; if (offset == stride-1) v += nn*stride; } } return m; } #define NO_INSTANTIATE_ARRAY_SORT_API(T, API) \ template <> API Array<T> \ Array<T>::sort (int, sortmode) const \ { \ return *this; \ } \ template <> API Array<T> \ Array<T>::sort (Array<octave_idx_type> &sidx, int, sortmode) const \ { \ sidx = Array<octave_idx_type> (); \ return *this; \ } \ template <> API sortmode \ Array<T>::issorted (sortmode) const \ { \ return UNSORTED; \ } \ API Array<T>::compare_fcn_type \ safe_comparator (sortmode, const Array<T>&, bool) \ { \ return nullptr; \ } \ template <> API Array<octave_idx_type> \ Array<T>::sort_rows_idx (sortmode) const \ { \ return Array<octave_idx_type> (); \ } \ template <> API sortmode \ Array<T>::is_sorted_rows (sortmode) const \ { \ return UNSORTED; \ } \ template <> API octave_idx_type \ Array<T>::lookup (T const &, sortmode) const \ { \ return 0; \ } \ template <> API Array<octave_idx_type> \ Array<T>::lookup (const Array<T>&, sortmode) const \ { \ return Array<octave_idx_type> (); \ } \ template <> API octave_idx_type \ Array<T>::nnz () const \ { \ return 0; \ } \ template <> API Array<octave_idx_type> \ Array<T>::find (octave_idx_type, bool) const \ { \ return Array<octave_idx_type> (); \ } \ template <> API Array<T> \ Array<T>::nth_element (const octave::idx_vector&, int) const { \ return Array<T> (); \ } #define NO_INSTANTIATE_ARRAY_SORT(T) NO_INSTANTIATE_ARRAY_SORT_API (T,) template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::diag (octave_idx_type k) const { dim_vector dv = dims (); octave_idx_type nd = dv.ndims (); Array<T, Alloc> d; if (nd > 2) (*current_liboctave_error_handler) ("Matrix must be 2-dimensional"); octave_idx_type nnr = dv(0); octave_idx_type nnc = dv(1); if (nnr == 0 && nnc == 0) ; // do nothing for empty matrix else if (nnr != 1 && nnc != 1) { // Extract diag from matrix if (k > 0) nnc -= k; else if (k < 0) nnr += k; if (nnr > 0 && nnc > 0) { octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; d.resize (dim_vector (ndiag, 1)); if (k > 0) { for (octave_idx_type i = 0; i < ndiag; i++) d.xelem (i) = elem (i, i+k); } else if (k < 0) { for (octave_idx_type i = 0; i < ndiag; i++) d.xelem (i) = elem (i-k, i); } else { for (octave_idx_type i = 0; i < ndiag; i++) d.xelem (i) = elem (i, i); } } else // Matlab returns [] 0x1 for out-of-range diagonal d.resize (dim_vector (0, 1)); } else { // Create diag matrix from vector octave_idx_type roff = 0; octave_idx_type coff = 0; if (k > 0) { roff = 0; coff = k; } else if (k < 0) { roff = -k; coff = 0; } if (nnr == 1) { octave_idx_type n = nnc + std::abs (k); d = Array<T, Alloc> (dim_vector (n, n), resize_fill_value ()); for (octave_idx_type i = 0; i < nnc; i++) d.xelem (i+roff, i+coff) = elem (0, i); } else { octave_idx_type n = nnr + std::abs (k); d = Array<T, Alloc> (dim_vector (n, n), resize_fill_value ()); for (octave_idx_type i = 0; i < nnr; i++) d.xelem (i+roff, i+coff) = elem (i, 0); } } return d; } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::diag (octave_idx_type m, octave_idx_type n) const { if (ndims () != 2 || (rows () != 1 && cols () != 1)) (*current_liboctave_error_handler) ("cat: invalid dimension"); Array<T, Alloc> retval (dim_vector (m, n), resize_fill_value ()); octave_idx_type nel = std::min (numel (), std::min (m, n)); for (octave_idx_type i = 0; i < nel; i++) retval.xelem (i, i) = xelem (i); return retval; } template <typename T, typename Alloc> Array<T, Alloc> Array<T, Alloc>::cat (int dim, octave_idx_type n, const Array<T, Alloc> *array_list) { // Default concatenation. bool (dim_vector::*concat_rule) (const dim_vector&, int) = &dim_vector::concat; if (dim == -1 || dim == -2) { concat_rule = &dim_vector::hvcat; dim = -dim - 1; } else if (dim < 0) (*current_liboctave_error_handler) ("cat: invalid dimension"); if (n == 1) return array_list[0]; else if (n == 0) return Array<T, Alloc> (); // Special case: // // cat (dim, [], ..., [], A, ...) // // with dim > 2, A not 0x0, and at least three arguments to // concatenate is equivalent to // // cat (dim, A, ...) // // Note that this check must be performed here because for full-on // braindead Matlab compatibility, we need to have things like // // cat (3, [], [], A) // // succeed, but to have things like // // cat (3, cat (3, [], []), A) // cat (3, zeros (0, 0, 2), A) // // fail. See also bug report #31615. octave_idx_type istart = 0; if (n > 2 && dim > 1) { for (octave_idx_type i = 0; i < n; i++) { dim_vector dv = array_list[i].dims (); if (dv.zero_by_zero ()) istart++; else break; } // Don't skip any initial arguments if they are all empty. if (istart >= n) istart = 0; } dim_vector dv = array_list[istart++].dims (); for (octave_idx_type i = istart; i < n; i++) if (! (dv.*concat_rule) (array_list[i].dims (), dim)) (*current_liboctave_error_handler) ("cat: dimension mismatch"); Array<T, Alloc> retval (dv); if (retval.isempty ()) return retval; int nidx = std::max (dv.ndims (), static_cast<octave_idx_type> (dim + 1)); Array<octave::idx_vector> idxa (dim_vector (nidx, 1), octave::idx_vector::colon); octave_idx_type l = 0; for (octave_idx_type i = 0; i < n; i++) { // NOTE: This takes some thinking, but no matter what the above rules // are, an empty array can always be skipped at this point, because // the result dimensions are already determined, and there is no way // an empty array may contribute a nonzero piece along the dimension // at this point, unless an empty array can be promoted to a non-empty // one (which makes no sense). I repeat, *no way*, think about it. if (array_list[i].isempty ()) continue; octave_quit (); octave_idx_type u; if (dim < array_list[i].ndims ()) u = l + array_list[i].dims ()(dim); else u = l + 1; idxa(dim) = octave::idx_vector (l, u); retval.assign (idxa, array_list[i]); l = u; } return retval; } template <typename T, typename Alloc> void Array<T, Alloc>::print_info (std::ostream& os, const std::string& prefix) const { os << prefix << "m_rep address: " << m_rep << '\n' << prefix << "m_rep->m_len: " << m_rep->m_len << '\n' << prefix << "m_rep->m_data: " << static_cast<void *> (m_rep->m_data) << '\n' << prefix << "m_rep->m_count: " << m_rep->m_count << '\n' << prefix << "m_slice_data: " << static_cast<void *> (m_slice_data) << '\n' << prefix << "m_slice_len: " << m_slice_len << '\n'; // 2D info: // // << pefix << "rows: " << rows () << "\n" // << prefix << "cols: " << cols () << "\n"; } template <typename T, typename Alloc> bool Array<T, Alloc>::optimize_dimensions (const dim_vector& dv) { bool retval = m_dimensions == dv; if (retval) m_dimensions = dv; return retval; } template <typename T, typename Alloc> void Array<T, Alloc>::instantiation_guard () { // This guards against accidental implicit instantiations. // Array<T, Alloc> instances should always be explicit and use INSTANTIATE_ARRAY. T::__xXxXx__ (); } #define INSTANTIATE_ARRAY(T, API) \ template <> API void \ Array<T>::instantiation_guard () { } \ \ template class API Array<T> // FIXME: is this used? template <typename T, typename Alloc> std::ostream& operator << (std::ostream& os, const Array<T, Alloc>& a) { dim_vector a_dims = a.dims (); int n_dims = a_dims.ndims (); os << n_dims << "-dimensional array"; if (n_dims) os << " (" << a_dims.str () << ')'; os <<"\n\n"; if (n_dims) { os << "data:"; Array<octave_idx_type> ra_idx (dim_vector (n_dims, 1), 0); // Number of times the first 2d-array is to be displayed. octave_idx_type m = 1; for (int i = 2; i < n_dims; i++) m *= a_dims(i); if (m == 1) { octave_idx_type rows = 0; octave_idx_type cols = 0; switch (n_dims) { case 2: rows = a_dims(0); cols = a_dims(1); for (octave_idx_type j = 0; j < rows; j++) { ra_idx(0) = j; for (octave_idx_type k = 0; k < cols; k++) { ra_idx(1) = k; os << ' ' << a.elem (ra_idx); } os << "\n"; } break; default: rows = a_dims(0); for (octave_idx_type k = 0; k < rows; k++) { ra_idx(0) = k; os << ' ' << a.elem (ra_idx); } break; } os << "\n"; } else { octave_idx_type rows = a_dims(0); octave_idx_type cols = a_dims(1); for (int i = 0; i < m; i++) { os << "\n(:,:,"; for (int j = 2; j < n_dims - 1; j++) os << ra_idx(j) + 1 << ','; os << ra_idx(n_dims - 1) + 1 << ") = \n"; for (octave_idx_type j = 0; j < rows; j++) { ra_idx(0) = j; for (octave_idx_type k = 0; k < cols; k++) { ra_idx(1) = k; os << ' ' << a.elem (ra_idx); } os << "\n"; } os << "\n"; if (i != m - 1) increment_index (ra_idx, a_dims, 2); } } } return os; }