view liboctave/array/Array-util.cc @ 31605:e88a07dec498 stable

maint: Use macros to begin/end C++ namespaces. * oct-conf-post-public.in.h: Define two macros (OCTAVE_BEGIN_NAMESPACE, OCTAVE_END_NAMESPACE) that can be used to start/end a namespace. * mk-opts.pl, build-env.h, build-env.in.cc, __betainc__.cc, __contourc__.cc, __dsearchn__.cc, __eigs__.cc, __expint__.cc, __ftp__.cc, __gammainc__.cc, __ichol__.cc, __ilu__.cc, __isprimelarge__.cc, __lin_interpn__.cc, __magick_read__.cc, __pchip_deriv__.cc, __qp__.cc, amd.cc, auto-shlib.cc, auto-shlib.h, balance.cc, base-text-renderer.cc, base-text-renderer.h, besselj.cc, bitfcns.cc, bsxfun.cc, c-file-ptr-stream.cc, c-file-ptr-stream.h, call-stack.cc, call-stack.h, ccolamd.cc, cellfun.cc, chol.cc, colamd.cc, colloc.cc, conv2.cc, daspk.cc, dasrt.cc, dassl.cc, data.cc, data.h, debug.cc, defaults.cc, defaults.h, defun-int.h, defun.cc, det.cc, dirfns.cc, display.cc, display.h, dlmread.cc, dmperm.cc, dot.cc, dynamic-ld.cc, dynamic-ld.h, eig.cc, ellipj.cc, environment.cc, environment.h, error.cc, error.h, errwarn.h, event-manager.cc, event-manager.h, event-queue.cc, event-queue.h, fcn-info.cc, fcn-info.h, fft.cc, fft2.cc, fftn.cc, file-io.cc, filter.cc, find.cc, ft-text-renderer.cc, ft-text-renderer.h, gcd.cc, getgrent.cc, getpwent.cc, getrusage.cc, givens.cc, gl-render.cc, gl-render.h, gl2ps-print.cc, gl2ps-print.h, graphics-toolkit.cc, graphics-toolkit.h, graphics.cc, graphics.in.h, gsvd.cc, gtk-manager.cc, gtk-manager.h, hash.cc, help.cc, help.h, hess.cc, hex2num.cc, hook-fcn.cc, hook-fcn.h, input.cc, input.h, interpreter-private.cc, interpreter-private.h, interpreter.cc, interpreter.h, inv.cc, jsondecode.cc, jsonencode.cc, kron.cc, latex-text-renderer.cc, latex-text-renderer.h, load-path.cc, load-path.h, load-save.cc, load-save.h, lookup.cc, ls-ascii-helper.cc, ls-ascii-helper.h, ls-oct-text.cc, ls-utils.cc, ls-utils.h, lsode.cc, lu.cc, mappers.cc, matrix_type.cc, max.cc, mex-private.h, mex.cc, mgorth.cc, nproc.cc, oct-fstrm.cc, oct-fstrm.h, oct-hdf5-types.cc, oct-hdf5-types.h, oct-hist.cc, oct-hist.h, oct-iostrm.cc, oct-iostrm.h, oct-opengl.h, oct-prcstrm.cc, oct-prcstrm.h, oct-procbuf.cc, oct-procbuf.h, oct-process.cc, oct-process.h, oct-stdstrm.h, oct-stream.cc, oct-stream.h, oct-strstrm.cc, oct-strstrm.h, oct-tex-lexer.in.ll, oct-tex-parser.yy, ordqz.cc, ordschur.cc, pager.cc, pager.h, pinv.cc, pow2.cc, pr-flt-fmt.cc, pr-output.cc, procstream.cc, procstream.h, psi.cc, qr.cc, quad.cc, quadcc.cc, qz.cc, rand.cc, rcond.cc, regexp.cc, schur.cc, settings.cc, settings.h, sighandlers.cc, sighandlers.h, sparse-xdiv.cc, sparse-xdiv.h, sparse-xpow.cc, sparse-xpow.h, sparse.cc, spparms.cc, sqrtm.cc, stack-frame.cc, stack-frame.h, stream-euler.cc, strfind.cc, strfns.cc, sub2ind.cc, svd.cc, sylvester.cc, symbfact.cc, syminfo.cc, syminfo.h, symrcm.cc, symrec.cc, symrec.h, symscope.cc, symscope.h, symtab.cc, symtab.h, syscalls.cc, sysdep.cc, sysdep.h, text-engine.cc, text-engine.h, text-renderer.cc, text-renderer.h, time.cc, toplev.cc, tril.cc, tsearch.cc, typecast.cc, url-handle-manager.cc, url-handle-manager.h, urlwrite.cc, utils.cc, utils.h, variables.cc, variables.h, xdiv.cc, xdiv.h, xnorm.cc, xnorm.h, xpow.cc, xpow.h, __delaunayn__.cc, __fltk_uigetfile__.cc, __glpk__.cc, __init_fltk__.cc, __init_gnuplot__.cc, __ode15__.cc, __voronoi__.cc, audiodevinfo.cc, audioread.cc, convhulln.cc, fftw.cc, gzip.cc, mk-build-env-features.sh, mk-builtins.pl, cdef-class.cc, cdef-class.h, cdef-fwd.h, cdef-manager.cc, 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.cc, ov-base.h, ov-bool-mat.cc, ov-builtin.h, ov-cell.cc, ov-class.cc, ov-class.h, ov-classdef.cc, ov-classdef.h, ov-complex.cc, ov-fcn-handle.cc, ov-fcn-handle.h, ov-fcn.h, ov-java.cc, ov-java.h, ov-mex-fcn.h, ov-null-mat.cc, ov-oncleanup.cc, ov-struct.cc, ov-typeinfo.cc, ov-typeinfo.h, ov-usr-fcn.cc, ov-usr-fcn.h, ov.cc, ov.h, octave.cc, octave.h, mk-ops.sh, op-b-b.cc, op-b-bm.cc, op-b-sbm.cc, op-bm-b.cc, op-bm-bm.cc, op-bm-sbm.cc, op-cdm-cdm.cc, op-cell.cc, op-chm.cc, op-class.cc, op-cm-cm.cc, op-cm-cs.cc, op-cm-m.cc, op-cm-s.cc, op-cm-scm.cc, op-cm-sm.cc, op-cs-cm.cc, op-cs-cs.cc, op-cs-m.cc, op-cs-s.cc, op-cs-scm.cc, op-cs-sm.cc, op-dm-dm.cc, op-dm-scm.cc, op-dm-sm.cc, op-dm-template.cc, op-dms-template.cc, op-fcdm-fcdm.cc, op-fcm-fcm.cc, op-fcm-fcs.cc, op-fcm-fm.cc, op-fcm-fs.cc, op-fcn.cc, op-fcs-fcm.cc, op-fcs-fcs.cc, op-fcs-fm.cc, op-fcs-fs.cc, op-fdm-fdm.cc, op-fm-fcm.cc, op-fm-fcs.cc, op-fm-fm.cc, op-fm-fs.cc, op-fs-fcm.cc, op-fs-fcs.cc, op-fs-fm.cc, op-fs-fs.cc, op-i16-i16.cc, op-i32-i32.cc, op-i64-i64.cc, op-i8-i8.cc, op-int-concat.cc, op-m-cm.cc, op-m-cs.cc, op-m-m.cc, op-m-s.cc, op-m-scm.cc, op-m-sm.cc, op-mi.cc, op-pm-pm.cc, op-pm-scm.cc, op-pm-sm.cc, op-pm-template.cc, op-range.cc, op-s-cm.cc, op-s-cs.cc, op-s-m.cc, op-s-s.cc, op-s-scm.cc, op-s-sm.cc, op-sbm-b.cc, op-sbm-bm.cc, op-sbm-sbm.cc, op-scm-cm.cc, op-scm-cs.cc, op-scm-m.cc, op-scm-s.cc, op-scm-scm.cc, op-scm-sm.cc, op-sm-cm.cc, op-sm-cs.cc, op-sm-m.cc, op-sm-s.cc, op-sm-scm.cc, op-sm-sm.cc, op-str-m.cc, op-str-s.cc, op-str-str.cc, op-struct.cc, op-ui16-ui16.cc, op-ui32-ui32.cc, op-ui64-ui64.cc, op-ui8-ui8.cc, ops.h, anon-fcn-validator.cc, 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.cc, 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.cc, pt-bp.h, pt-cbinop.cc, pt-cbinop.h, pt-cell.cc, pt-cell.h, pt-check.cc, pt-check.h, pt-classdef.cc, pt-classdef.h, pt-cmd.h, pt-colon.cc, pt-colon.h, pt-const.cc, 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.cc, 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-vm-eval.cc, pt-walk.cc, pt-walk.h, pt.cc, pt.h, token.cc, token.h, Range.cc, Range.h, idx-vector.cc, idx-vector.h, range-fwd.h, CollocWt.cc, CollocWt.h, aepbalance.cc, aepbalance.h, chol.cc, chol.h, gepbalance.cc, gepbalance.h, gsvd.cc, gsvd.h, hess.cc, hess.h, lo-mappers.cc, lo-mappers.h, lo-specfun.cc, lo-specfun.h, lu.cc, lu.h, oct-convn.cc, oct-convn.h, oct-fftw.cc, oct-fftw.h, oct-norm.cc, oct-norm.h, oct-rand.cc, oct-rand.h, oct-spparms.cc, oct-spparms.h, qr.cc, qr.h, qrp.cc, qrp.h, randgamma.cc, randgamma.h, randmtzig.cc, randmtzig.h, randpoisson.cc, randpoisson.h, schur.cc, 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.cc, action-container.h, base-list.h, cmd-edit.cc, cmd-edit.h, cmd-hist.cc, cmd-hist.h, f77-fcn.h, file-info.cc, file-info.h, lo-array-errwarn.cc, lo-array-errwarn.h, lo-hash.cc, lo-hash.h, lo-ieee.h, lo-regexp.cc, lo-regexp.h, lo-utils.cc, lo-utils.h, oct-base64.cc, oct-base64.h, oct-glob.cc, oct-glob.h, oct-inttypes.h, oct-mutex.cc, oct-mutex.h, oct-refcount.h, oct-shlib.cc, oct-shlib.h, oct-sparse.cc, oct-sparse.h, oct-string.h, octave-preserve-stream-state.h, pathsearch.cc, pathsearch.h, quit.cc, quit.h, unwind-prot.cc, unwind-prot.h, url-transfer.cc, url-transfer.h : Use new macros to begin/end C++ namespaces.
author Rik <rik@octave.org>
date Thu, 01 Dec 2022 14:23:45 -0800
parents 972959edc3ff
children 597f3ee61a48
line wrap: on
line source

////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2003-2022 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/>.
//
////////////////////////////////////////////////////////////////////////

#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif

#include "Array-util.h"
#include "lo-error.h"
#include "oct-locbuf.h"

bool
index_in_bounds (const Array<octave_idx_type>& ra_idx,
                 const dim_vector& dimensions)
{
  bool retval = true;

  int n = ra_idx.numel ();

  if (n == dimensions.ndims ())
    {
      for (int i = 0; i < n; i++)
        {
          if (ra_idx(i) < 0 || ra_idx(i) >= dimensions(i))
            {
              retval = false;
              break;
            }
        }
    }
  else
    retval = false;

  return retval;
}

void
increment_index (Array<octave_idx_type>& ra_idx, const dim_vector& dimensions,
                 int start_dimension)
{
  ra_idx(start_dimension)++;

  int n = ra_idx.numel () - 1;
  int nda = dimensions.ndims ();

  for (int i = start_dimension; i < n; i++)
    {
      if (ra_idx(i) < (i < nda ? dimensions(i) : 1))
        break;
      else
        {
          ra_idx(i) = 0;
          ra_idx(i+1)++;
        }
    }
}

octave_idx_type
get_scalar_idx (Array<octave_idx_type>& idx, dim_vector& dims)
{
  octave_idx_type retval (-1);

  int n = idx.numel ();

  if (n > 0)
    {
      retval = idx(--n);

      while (--n >= 0)
        {
          retval *= dims(n);

          retval += idx(n);
        }
    }
  return retval;
}

octave_idx_type
num_ones (const Array<octave_idx_type>& ra_idx)
{
  octave_idx_type retval = 0;

  for (octave_idx_type i = 0; i < ra_idx.numel (); i++)
    {
      if (ra_idx(i) == 1)
        retval++;
    }

  return retval;
}

bool
is_scalar (const dim_vector& dim)
{
  bool retval = true;

  int n = dim.ndims ();

  if (n == 0)
    retval = false;
  else
    {
      for (int i = 0; i < n; i++)
        {
          if (dim(i) != 1)
            {
              retval = false;

              break;
            }
        }
    }
  return retval;
}

bool
isvector (const dim_vector& dim)
{
  int m = 0;
  int n = dim.ndims ();

  if (n == 0)
    m = 2;
  else
    {
      for (int i = 0; i < n; i++)
        if (dim(i) > 1)
          m++;
        else if (dim(i) < 1)
          m += 2;
    }

  return (m < 2);
}

bool
any_ones (const Array<octave_idx_type>& arr)
{
  bool retval = false;

  for (octave_idx_type i = 0; i < arr.numel (); i++)
    {
      if (arr (i) == 1)
        {
          retval = true;

          break;
        }
    }
  return retval;
}

octave_idx_type
compute_index (octave_idx_type n, const dim_vector& dims)
{
  if (n < 0)
    octave::err_invalid_index (n, 1, 1);
  if (n >= dims.numel ())
    octave::err_index_out_of_range (1, 1, n+1, dims.numel (), dims);

  return n;
}

octave_idx_type
compute_index (octave_idx_type i, octave_idx_type j, const dim_vector& dims)
{
  if (i < 0)
    octave::err_invalid_index (i, 2, 1);
  if (j < 0)
    octave::err_invalid_index (j, 2, 2);
  if (i >= dims(0))
    octave::err_index_out_of_range (2, 1, i+1, dims(0), dims);
  if (j >= dims.numel (1))
    octave::err_index_out_of_range (2, 2, j+1, dims.numel (1), dims);

  return j*dims(0) + i;
}

octave_idx_type
compute_index (octave_idx_type i, octave_idx_type j, octave_idx_type k,
               const dim_vector& dims)
{
  if (i < 0)
    octave::err_invalid_index (i, 3, 1);
  if (j < 0)
    octave::err_invalid_index (j, 3, 2);
  if (k < 0)
    octave::err_invalid_index (k, 3, 3);
  if (i >= dims(0))
    octave::err_index_out_of_range (3, 1, i+1, dims(0), dims);
  if (j >= dims(1))
    octave::err_index_out_of_range (3, 2, j+1, dims(1), dims);
  if (k >= dims.numel (2))
    octave::err_index_out_of_range (3, 3, k+1, dims.numel (2), dims);

  return (k*dims(1) + j)*dims(0) + i;
}

octave_idx_type
compute_index (const Array<octave_idx_type>& ra_idx, const dim_vector& dims)
{
  int nd = ra_idx.numel ();
  const dim_vector dv = dims.redim (nd);
  for (int d = 0; d < nd; d++)
    {
      if (ra_idx(d) < 0)
        octave::err_invalid_index (ra_idx(d), nd, d+1);
      if (ra_idx(d) >= dv(d))
        octave::err_index_out_of_range (nd, d+1, ra_idx(d)+1, dv(d), dims);
    }

  return dv.compute_index (ra_idx.data ());
}

Array<octave_idx_type>
conv_to_int_array (const Array<octave::idx_vector>& a)
{
  Array<octave_idx_type> retval (a.dims ());

  for (octave_idx_type i = 0; i < a.numel (); i++)
    retval(i) = a(i).elem (0);

  return retval;
}

Array<octave::idx_vector>
conv_to_array (const octave::idx_vector *tmp, const octave_idx_type len)
{
  Array<octave::idx_vector> retval (dim_vector (len, 1));

  for (octave_idx_type i = 0; i < len; i++)
    retval(i) = tmp[i];

  return retval;
}

dim_vector
freeze (Array<octave::idx_vector>& ra_idx, const dim_vector& dimensions, int resize_ok)
{
  dim_vector retval;

  int n = ra_idx.numel ();

  assert (n == dimensions.ndims ());

  retval.resize (n);

  static const char *tag[3] = { "row", "column", nullptr };

  for (int i = 0; i < n; i++)
    retval(i) = ra_idx(i).freeze (dimensions(i), tag[i < 2 ? i : 2],
                                  resize_ok);

  return retval;
}

bool
vector_equivalent (const dim_vector& dv)
{
  int n = dv.ndims ();

  bool found_first = false;

  for (int i = 0; i < n; i++)
    {
      if (dv(i) != 1)
        {
          if (! found_first)
            found_first = true;
          else
            return false;
        }
    }

  return true;
}

bool
all_ok (const Array<octave::idx_vector>& ra_idx)
{
  bool retval = true;

  octave_idx_type n = ra_idx.numel ();

  for (octave_idx_type i = 0; i < n; i++)
    {
      if (! ra_idx(i))
        {
          retval = false;
          break;
        }
    }

  return retval;
}

bool
any_orig_empty (const Array<octave::idx_vector>& ra_idx)
{
  bool retval = false;

  octave_idx_type n = ra_idx.numel ();

  for (octave_idx_type i = 0; i < n; i++)
    {
      if (ra_idx(i).orig_empty ())
        {
          retval = true;
          break;
        }
    }

  return retval;
}

bool
all_colon_equiv (const Array<octave::idx_vector>& ra_idx,
                 const dim_vector& frozen_lengths)
{
  bool retval = true;

  octave_idx_type idx_n = ra_idx.numel ();

  int n = frozen_lengths.ndims ();

  assert (idx_n == n);

  for (octave_idx_type i = 0; i < n; i++)
    {
      if (! ra_idx(i).is_colon_equiv (frozen_lengths(i)))
        {
          retval = false;
          break;
        }
    }

  return retval;
}

bool
all_ones (const Array<octave_idx_type>& arr)
{
  bool retval = true;

  for (octave_idx_type i = 0; i < arr.numel (); i++)
    {
      if (arr(i) != 1)
        {
          retval = false;
          break;
        }
    }

  return retval;
}

Array<octave_idx_type>
get_elt_idx (const Array<octave::idx_vector>& ra_idx,
             const Array<octave_idx_type>& result_idx)
{
  octave_idx_type n = ra_idx.numel ();

  Array<octave_idx_type> retval (dim_vector (n, 1));

  for (octave_idx_type i = 0; i < n; i++)
    retval(i) = ra_idx(i).elem (result_idx(i));

  return retval;
}

Array<octave_idx_type>
get_ra_idx (octave_idx_type idx, const dim_vector& dims)
{
  Array<octave_idx_type> retval;

  int n_dims = dims.ndims ();

  retval.resize (dim_vector (n_dims, 1));

  for (int i = 0; i < n_dims; i++)
    retval(i) = 0;

  assert (idx > 0 || idx < dims.numel ());

  for (octave_idx_type i = 0; i < idx; i++)
    increment_index (retval, dims);

  // FIXME: the solution using increment_index is not efficient.

#if 0
  octave_idx_type var = 1;
  for (int i = 0; i < n_dims; i++)
    {
      std::cout << "idx: " << idx << ", var: " << var
                << ", dims(" << i << "): " << dims(i) <<"\n";
      retval(i) = ((int)floor(((idx) / (double)var))) % dims(i);
      idx -= var * retval(i);
      var = dims(i);
    }
#endif

  return retval;
}

dim_vector
zero_dims_inquire (const Array<octave::idx_vector>& ia, const dim_vector& rhdv)
{
  int ial = ia.numel ();
  int rhdvl = rhdv.ndims ();
  dim_vector rdv = dim_vector::alloc (ial);
  bool *scalar = new bool [ial];
  bool *colon = new bool [ial];
  // Mark scalars and colons, count non-scalar indices.
  int nonsc = 0;
  bool all_colons = true;
  for (int i = 0; i < ial; i++)
    {
      // FIXME: should we check for length() instead?
      scalar[i] = ia(i).is_scalar ();
      colon[i] = ia(i).is_colon ();
      if (! scalar[i]) nonsc++;
      if (! colon[i]) rdv(i) = ia(i).extent (0);
      all_colons = all_colons && colon[i];
    }

  // If the number of nonscalar indices matches the dimensionality of
  // RHS, we try an exact match, inquiring even singleton dimensions.
  if (all_colons)
    {
      rdv = rhdv;
      rdv.resize (ial, 1);
    }
  else if (nonsc == rhdvl)
    {
      for (int i = 0, j = 0; i < ial; i++)
        {
          if (scalar[i]) continue;
          if (colon[i])
            rdv(i) = rhdv(j);
          j++;
        }
    }
  else
    {
      dim_vector rhdv0 = rhdv;
      rhdv0.chop_all_singletons ();
      int rhdv0l = rhdv0.ndims ();
      for (int i = 0, j = 0; i < ial; i++)
        {
          if (scalar[i]) continue;
          if (colon[i])
            rdv(i) = (j < rhdv0l) ? rhdv0(j++) : 1;
        }
    }

  delete [] scalar;
  delete [] colon;

  return rdv;
}

dim_vector
zero_dims_inquire (const octave::idx_vector& i, const octave::idx_vector& j,
                   const dim_vector& rhdv)
{
  bool icol = i.is_colon ();
  bool jcol = j.is_colon ();
  dim_vector rdv;
  if (icol && jcol && rhdv.ndims () == 2)
    {
      rdv(0) = rhdv(0);
      rdv(1) = rhdv(1);
    }
  else if (rhdv.ndims () == 2
           && ! i.is_scalar () && ! j.is_scalar ())
    {
      rdv(0) = (icol ? rhdv(0) : i.extent (0));
      rdv(1) = (jcol ? rhdv(1) : j.extent (0));
    }
  else
    {
      dim_vector rhdv0 = rhdv;
      rhdv0.chop_all_singletons ();
      int k = 0;
      rdv(0) = i.extent (0);
      if (icol)
        rdv(0) = rhdv0(k++);
      else if (! i.is_scalar ())
        k++;
      rdv(1) = j.extent (0);
      if (jcol)
        rdv(1) = rhdv0(k++);
      else if (! j.is_scalar ())
        k++;
    }

  return rdv;
}

octave::idx_vector
sub2ind (const dim_vector& dv, const Array<octave::idx_vector>& idxa)
{
  octave::idx_vector retval;
  octave_idx_type len = idxa.numel ();

  if (len == 0)
    (*current_liboctave_error_handler) ("sub2ind: needs at least 2 indices");

  const dim_vector dvx = dv.redim (len);
  bool all_ranges = true;
  octave_idx_type clen = -1;

  for (octave_idx_type i = 0; i < len; i++)
    {
      try
        {
          octave::idx_vector idx = idxa(i);
          octave_idx_type n = dvx(i);

          all_ranges = all_ranges && idx.is_range ();
          if (clen < 0)
            clen = idx.length (n);
          else if (clen != idx.length (n))
            (*current_liboctave_error_handler)
              ("sub2ind: lengths of indices must match");

          if (idx.extent (n) > n)
            octave::err_index_out_of_range (len, i+1, idx.extent (n), n, dv);
        }
      catch (octave::index_exception& ie)
        {
          ie.set_pos_if_unset (len, i+1);
          ie.set_var ();
          std::string msg = ie.message ();
          (*current_liboctave_error_with_id_handler)
            (ie.err_id (), "%s", msg.c_str ());
        }
    }
  // idxa known to be valid.
  // Shouldn't need to catch index_exception below here.

  if (len == 1)
    retval = idxa(0);
  else if (clen == 1)
    {
      // All scalars case - the result is a scalar.
      octave_idx_type idx = idxa(len-1)(0);
      for (octave_idx_type i = len - 2; i >= 0; i--)
        idx = dvx(i) * idx + idxa(i)(0);
      retval = octave::idx_vector (idx);
    }
  else if (all_ranges && clen != 0)
    {
      // All ranges case - the result is a range.
      octave_idx_type start = 0;
      octave_idx_type step = 0;
      for (octave_idx_type i = len - 1; i >= 0; i--)
        {
          octave_idx_type xstart = idxa(i)(0);
          octave_idx_type xstep = idxa(i)(1) - xstart;
          start = dvx(i) * start + xstart;
          step = dvx(i) * step + xstep;
        }
      retval = octave::idx_vector::make_range (start, step, clen);
    }
  else
    {
      Array<octave_idx_type> idx (idxa(0).orig_dimensions ());

      for (octave_idx_type i = len - 1; i >= 0; i--)
        {
          // Initialized inside the loop so that each call to
          // idx_vector::loop operates from the beginning of IDX_VEC.

          octave_idx_type *idx_vec = idx.fortran_vec ();

          if (i < len - 1)
            {
              octave_idx_type n = dvx(i);

              idxa(i).loop (clen, [=, &idx_vec] (octave_idx_type k) {
                (*idx_vec++ *= n) += k;
              });
            }
          else
            idxa(i).copy_data (idx_vec);
        }

      retval = octave::idx_vector (idx);
    }

  return retval;
}

Array<octave::idx_vector>
ind2sub (const dim_vector& dv, const octave::idx_vector& idx)
{
  octave_idx_type len = idx.length (0);
  octave_idx_type n = dv.ndims ();
  Array<octave::idx_vector> retval (dim_vector (n, 1));
  octave_idx_type numel = dv.numel ();

  if (idx.is_scalar ())
    {
      octave_idx_type k = idx(0);
      for (octave_idx_type j = 0; j < n-1; j++)
        {
          retval(j) = k % dv(j);
          k /= dv(j);
        }

      retval(n-1) = idx(0) < numel ? k % dv(n-1) : k;
    }
  else
    {
      OCTAVE_LOCAL_BUFFER (Array<octave_idx_type>, rdata, n);

      dim_vector odv = idx.orig_dimensions ();
      for (octave_idx_type j = 0; j < n; j++)
        rdata[j] = Array<octave_idx_type> (odv);

      for (octave_idx_type i = 0; i < len; i++)
        {
          octave_idx_type k = idx(i);
          for (octave_idx_type j = 0; j < n-1; j++)
            {
              rdata[j](i) = k % dv(j);
              k /= dv(j);
            }

          rdata[n-1](i) = idx(i) < numel ? k % dv(n-1) : k;
        }

      for (octave_idx_type j = 0; j < n; j++)
        retval(j) = rdata[j];
    }

  return retval;
}

int
permute_vector_compare (const void *a, const void *b)
{
  const permute_vector *pva = static_cast<const permute_vector *> (a);
  const permute_vector *pvb = static_cast<const permute_vector *> (b);

  return pva->pidx > pvb->pidx;
}