Mercurial > hg > octave-nkf

comparison liboctave/array/Array.cc @ 15271:648dabbb4c6b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
build: Refactor liboctave into multiple subdirectories. Move libcruft into liboctave. * array/Array-C.cc, array/Array-b.cc, array/Array-ch.cc, array/Array-d.cc, array/Array-f.cc, array/Array-fC.cc, array/Array-i.cc, array/Array-idx-vec.cc, array/Array-s.cc, array/Array-str.cc, array/Array-util.cc, array/Array-util.h, array/Array-voidp.cc, array/Array.cc, array/Array.h, array/Array2.h, array/Array3.h, array/ArrayN.h, array/CColVector.cc, array/CColVector.h, array/CDiagMatrix.cc, array/CDiagMatrix.h, array/CMatrix.cc, array/CMatrix.h, array/CNDArray.cc, array/CNDArray.h, array/CRowVector.cc, array/CRowVector.h, array/CSparse.cc, array/CSparse.h, array/DiagArray2.cc, array/DiagArray2.h, array/MArray-C.cc, array/MArray-d.cc, array/MArray-decl.h, array/MArray-defs.h, array/MArray-f.cc, array/MArray-fC.cc, array/MArray-i.cc, array/MArray-s.cc, array/MArray.cc, array/MArray.h, array/MArray2.h, array/MArrayN.h, array/MDiagArray2.cc, array/MDiagArray2.h, array/MSparse-C.cc, array/MSparse-d.cc, array/MSparse-defs.h, array/MSparse.cc, array/MSparse.h, array/Matrix.h, array/MatrixType.cc, array/MatrixType.h, array/PermMatrix.cc, array/PermMatrix.h, array/Range.cc, array/Range.h, array/Sparse-C.cc, array/Sparse-b.cc, array/Sparse-d.cc, array/Sparse.cc, array/Sparse.h, array/boolMatrix.cc, array/boolMatrix.h, array/boolNDArray.cc, array/boolNDArray.h, array/boolSparse.cc, array/boolSparse.h, array/chMatrix.cc, array/chMatrix.h, array/chNDArray.cc, array/chNDArray.h, array/dColVector.cc, array/dColVector.h, array/dDiagMatrix.cc, array/dDiagMatrix.h, array/dMatrix.cc, array/dMatrix.h, array/dNDArray.cc, array/dNDArray.h, array/dRowVector.cc, array/dRowVector.h, array/dSparse.cc, array/dSparse.h, array/dim-vector.cc, array/dim-vector.h, array/fCColVector.cc, array/fCColVector.h, array/fCDiagMatrix.cc, array/fCDiagMatrix.h, array/fCMatrix.cc, array/fCMatrix.h, array/fCNDArray.cc, array/fCNDArray.h, array/fCRowVector.cc, array/fCRowVector.h, array/fColVector.cc, array/fColVector.h, array/fDiagMatrix.cc, array/fDiagMatrix.h, array/fMatrix.cc, array/fMatrix.h, array/fNDArray.cc, array/fNDArray.h, array/fRowVector.cc, array/fRowVector.h, array/idx-vector.cc, array/idx-vector.h, array/int16NDArray.cc, array/int16NDArray.h, array/int32NDArray.cc, array/int32NDArray.h, array/int64NDArray.cc, array/int64NDArray.h, array/int8NDArray.cc, array/int8NDArray.h, array/intNDArray.cc, array/intNDArray.h, array/module.mk, array/uint16NDArray.cc, array/uint16NDArray.h, array/uint32NDArray.cc, array/uint32NDArray.h, array/uint64NDArray.cc, array/uint64NDArray.h, array/uint8NDArray.cc, array/uint8NDArray.h: Moved from liboctave dir to array subdirectory. * cruft/Makefile.am, cruft/amos/README, cruft/amos/cacai.f, cruft/amos/cacon.f, cruft/amos/cairy.f, cruft/amos/casyi.f, cruft/amos/cbesh.f, cruft/amos/cbesi.f, cruft/amos/cbesj.f, cruft/amos/cbesk.f, cruft/amos/cbesy.f, cruft/amos/cbinu.f, cruft/amos/cbiry.f, cruft/amos/cbknu.f, cruft/amos/cbuni.f, cruft/amos/cbunk.f, cruft/amos/ckscl.f, cruft/amos/cmlri.f, cruft/amos/crati.f, cruft/amos/cs1s2.f, cruft/amos/cseri.f, cruft/amos/cshch.f, cruft/amos/cuchk.f, cruft/amos/cunhj.f, cruft/amos/cuni1.f, cruft/amos/cuni2.f, cruft/amos/cunik.f, cruft/amos/cunk1.f, cruft/amos/cunk2.f, cruft/amos/cuoik.f, cruft/amos/cwrsk.f, cruft/amos/dgamln.f, cruft/amos/gamln.f, cruft/amos/module.mk, cruft/amos/xzabs.f, cruft/amos/xzexp.f, cruft/amos/xzlog.f, cruft/amos/xzsqrt.f, cruft/amos/zacai.f, cruft/amos/zacon.f, cruft/amos/zairy.f, cruft/amos/zasyi.f, cruft/amos/zbesh.f, cruft/amos/zbesi.f, cruft/amos/zbesj.f, cruft/amos/zbesk.f, cruft/amos/zbesy.f, cruft/amos/zbinu.f, cruft/amos/zbiry.f, cruft/amos/zbknu.f, cruft/amos/zbuni.f, cruft/amos/zbunk.f, cruft/amos/zdiv.f, cruft/amos/zkscl.f, cruft/amos/zmlri.f, cruft/amos/zmlt.f, cruft/amos/zrati.f, cruft/amos/zs1s2.f, cruft/amos/zseri.f, cruft/amos/zshch.f, cruft/amos/zuchk.f, cruft/amos/zunhj.f, cruft/amos/zuni1.f, cruft/amos/zuni2.f, cruft/amos/zunik.f, cruft/amos/zunk1.f, cruft/amos/zunk2.f, cruft/amos/zuoik.f, cruft/amos/zwrsk.f, cruft/blas-xtra/cconv2.f, cruft/blas-xtra/cdotc3.f, cruft/blas-xtra/cmatm3.f, cruft/blas-xtra/csconv2.f, cruft/blas-xtra/dconv2.f, cruft/blas-xtra/ddot3.f, cruft/blas-xtra/dmatm3.f, cruft/blas-xtra/module.mk, cruft/blas-xtra/sconv2.f, cruft/blas-xtra/sdot3.f, cruft/blas-xtra/smatm3.f, cruft/blas-xtra/xcdotc.f, cruft/blas-xtra/xcdotu.f, cruft/blas-xtra/xddot.f, cruft/blas-xtra/xdnrm2.f, cruft/blas-xtra/xdznrm2.f, cruft/blas-xtra/xerbla.f, cruft/blas-xtra/xscnrm2.f, cruft/blas-xtra/xsdot.f, cruft/blas-xtra/xsnrm2.f, cruft/blas-xtra/xzdotc.f, cruft/blas-xtra/xzdotu.f, cruft/blas-xtra/zconv2.f, cruft/blas-xtra/zdconv2.f, cruft/blas-xtra/zdotc3.f, cruft/blas-xtra/zmatm3.f, cruft/daspk/datv.f, cruft/daspk/dcnst0.f, cruft/daspk/dcnstr.f, cruft/daspk/ddasic.f, cruft/daspk/ddasid.f, cruft/daspk/ddasik.f, cruft/daspk/ddaspk.f, cruft/daspk/ddstp.f, cruft/daspk/ddwnrm.f, cruft/daspk/dfnrmd.f, cruft/daspk/dfnrmk.f, cruft/daspk/dhels.f, cruft/daspk/dheqr.f, cruft/daspk/dinvwt.f, cruft/daspk/dlinsd.f, cruft/daspk/dlinsk.f, cruft/daspk/dmatd.f, cruft/daspk/dnedd.f, cruft/daspk/dnedk.f, cruft/daspk/dnsd.f, cruft/daspk/dnsid.f, cruft/daspk/dnsik.f, cruft/daspk/dnsk.f, cruft/daspk/dorth.f, cruft/daspk/dslvd.f, cruft/daspk/dslvk.f, cruft/daspk/dspigm.f, cruft/daspk/dyypnw.f, cruft/daspk/module.mk, cruft/dasrt/ddasrt.f, cruft/dasrt/drchek.f, cruft/dasrt/droots.f, cruft/dasrt/module.mk, cruft/dassl/ddaini.f, cruft/dassl/ddajac.f, cruft/dassl/ddanrm.f, cruft/dassl/ddaslv.f, cruft/dassl/ddassl.f, cruft/dassl/ddastp.f, cruft/dassl/ddatrp.f, cruft/dassl/ddawts.f, cruft/dassl/module.mk, cruft/fftpack/cfftb.f, cruft/fftpack/cfftb1.f, cruft/fftpack/cfftf.f, cruft/fftpack/cfftf1.f, cruft/fftpack/cffti.f, cruft/fftpack/cffti1.f, cruft/fftpack/fftpack.doc, cruft/fftpack/module.mk, cruft/fftpack/passb.f, cruft/fftpack/passb2.f, cruft/fftpack/passb3.f, cruft/fftpack/passb4.f, cruft/fftpack/passb5.f, cruft/fftpack/passf.f, cruft/fftpack/passf2.f, cruft/fftpack/passf3.f, cruft/fftpack/passf4.f, cruft/fftpack/passf5.f, cruft/fftpack/zfftb.f, cruft/fftpack/zfftb1.f, cruft/fftpack/zfftf.f, cruft/fftpack/zfftf1.f, cruft/fftpack/zffti.f, cruft/fftpack/zffti1.f, cruft/fftpack/zpassb.f, cruft/fftpack/zpassb2.f, cruft/fftpack/zpassb3.f, cruft/fftpack/zpassb4.f, cruft/fftpack/zpassb5.f, cruft/fftpack/zpassf.f, cruft/fftpack/zpassf2.f, cruft/fftpack/zpassf3.f, cruft/fftpack/zpassf4.f, cruft/fftpack/zpassf5.f, cruft/lapack-xtra/crsf2csf.f, cruft/lapack-xtra/module.mk, cruft/lapack-xtra/xclange.f, cruft/lapack-xtra/xdlamch.f, cruft/lapack-xtra/xdlange.f, cruft/lapack-xtra/xilaenv.f, cruft/lapack-xtra/xslamch.f, cruft/lapack-xtra/xslange.f, cruft/lapack-xtra/xzlange.f, cruft/lapack-xtra/zrsf2csf.f, cruft/link-deps.mk, cruft/misc/blaswrap.c, cruft/misc/cquit.c, cruft/misc/d1mach-tst.for, cruft/misc/d1mach.f, cruft/misc/f77-extern.cc, cruft/misc/f77-fcn.c, cruft/misc/f77-fcn.h, cruft/misc/i1mach.f, cruft/misc/lo-error.c, cruft/misc/lo-error.h, cruft/misc/module.mk, cruft/misc/quit.cc, cruft/misc/quit.h, cruft/misc/r1mach.f, cruft/mkf77def.in, cruft/odepack/cfode.f, cruft/odepack/dlsode.f, cruft/odepack/ewset.f, cruft/odepack/intdy.f, cruft/odepack/module.mk, cruft/odepack/prepj.f, cruft/odepack/scfode.f, cruft/odepack/sewset.f, cruft/odepack/sintdy.f, cruft/odepack/slsode.f, cruft/odepack/solsy.f, cruft/odepack/sprepj.f, cruft/odepack/ssolsy.f, cruft/odepack/sstode.f, cruft/odepack/stode.f, cruft/odepack/svnorm.f, cruft/odepack/vnorm.f, cruft/ordered-qz/README, cruft/ordered-qz/dsubsp.f, cruft/ordered-qz/exchqz.f, cruft/ordered-qz/module.mk, cruft/ordered-qz/sexchqz.f, cruft/ordered-qz/ssubsp.f, cruft/quadpack/dqagi.f, cruft/quadpack/dqagie.f, cruft/quadpack/dqagp.f, cruft/quadpack/dqagpe.f, cruft/quadpack/dqelg.f, cruft/quadpack/dqk15i.f, cruft/quadpack/dqk21.f, cruft/quadpack/dqpsrt.f, cruft/quadpack/module.mk, cruft/quadpack/qagi.f, cruft/quadpack/qagie.f, cruft/quadpack/qagp.f, cruft/quadpack/qagpe.f, cruft/quadpack/qelg.f, cruft/quadpack/qk15i.f, cruft/quadpack/qk21.f, cruft/quadpack/qpsrt.f, cruft/quadpack/xerror.f, cruft/ranlib/Basegen.doc, cruft/ranlib/HOWTOGET, cruft/ranlib/README, cruft/ranlib/advnst.f, cruft/ranlib/genbet.f, cruft/ranlib/genchi.f, cruft/ranlib/genexp.f, cruft/ranlib/genf.f, cruft/ranlib/gengam.f, cruft/ranlib/genmn.f, cruft/ranlib/genmul.f, cruft/ranlib/gennch.f, cruft/ranlib/gennf.f, cruft/ranlib/gennor.f, cruft/ranlib/genprm.f, cruft/ranlib/genunf.f, cruft/ranlib/getcgn.f, cruft/ranlib/getsd.f, cruft/ranlib/ignbin.f, cruft/ranlib/ignlgi.f, cruft/ranlib/ignnbn.f, cruft/ranlib/ignpoi.f, cruft/ranlib/ignuin.f, cruft/ranlib/initgn.f, cruft/ranlib/inrgcm.f, cruft/ranlib/lennob.f, cruft/ranlib/mltmod.f, cruft/ranlib/module.mk, cruft/ranlib/phrtsd.f, cruft/ranlib/qrgnin.f, cruft/ranlib/randlib.chs, cruft/ranlib/randlib.fdoc, cruft/ranlib/ranf.f, cruft/ranlib/setall.f, cruft/ranlib/setant.f, cruft/ranlib/setgmn.f, cruft/ranlib/setsd.f, cruft/ranlib/sexpo.f, cruft/ranlib/sgamma.f, cruft/ranlib/snorm.f, cruft/ranlib/tstbot.for, cruft/ranlib/tstgmn.for, cruft/ranlib/tstmid.for, cruft/ranlib/wrap.f, cruft/slatec-err/fdump.f, cruft/slatec-err/ixsav.f, cruft/slatec-err/j4save.f, cruft/slatec-err/module.mk, cruft/slatec-err/xerclr.f, cruft/slatec-err/xercnt.f, cruft/slatec-err/xerhlt.f, cruft/slatec-err/xermsg.f, cruft/slatec-err/xerprn.f, cruft/slatec-err/xerrwd.f, cruft/slatec-err/xersve.f, cruft/slatec-err/xgetf.f, cruft/slatec-err/xgetua.f, cruft/slatec-err/xsetf.f, cruft/slatec-err/xsetua.f, cruft/slatec-fn/acosh.f, cruft/slatec-fn/albeta.f, cruft/slatec-fn/algams.f, cruft/slatec-fn/alngam.f, cruft/slatec-fn/alnrel.f, cruft/slatec-fn/asinh.f, cruft/slatec-fn/atanh.f, cruft/slatec-fn/betai.f, cruft/slatec-fn/csevl.f, cruft/slatec-fn/d9gmit.f, cruft/slatec-fn/d9lgic.f, cruft/slatec-fn/d9lgit.f, cruft/slatec-fn/d9lgmc.f, cruft/slatec-fn/dacosh.f, cruft/slatec-fn/dasinh.f, cruft/slatec-fn/datanh.f, cruft/slatec-fn/dbetai.f, cruft/slatec-fn/dcsevl.f, cruft/slatec-fn/derf.f, cruft/slatec-fn/derfc.in.f, cruft/slatec-fn/dgami.f, cruft/slatec-fn/dgamit.f, cruft/slatec-fn/dgamlm.f, cruft/slatec-fn/dgamma.f, cruft/slatec-fn/dgamr.f, cruft/slatec-fn/dlbeta.f, cruft/slatec-fn/dlgams.f, cruft/slatec-fn/dlngam.f, cruft/slatec-fn/dlnrel.f, cruft/slatec-fn/dpchim.f, cruft/slatec-fn/dpchst.f, cruft/slatec-fn/erf.f, cruft/slatec-fn/erfc.in.f, cruft/slatec-fn/gami.f, cruft/slatec-fn/gamit.f, cruft/slatec-fn/gamlim.f, cruft/slatec-fn/gamma.f, cruft/slatec-fn/gamr.f, cruft/slatec-fn/initds.f, cruft/slatec-fn/inits.f, cruft/slatec-fn/module.mk, cruft/slatec-fn/pchim.f, cruft/slatec-fn/pchst.f, cruft/slatec-fn/r9gmit.f, cruft/slatec-fn/r9lgic.f, cruft/slatec-fn/r9lgit.f, cruft/slatec-fn/r9lgmc.f, cruft/slatec-fn/xacosh.f, cruft/slatec-fn/xasinh.f, cruft/slatec-fn/xatanh.f, cruft/slatec-fn/xbetai.f, cruft/slatec-fn/xdacosh.f, cruft/slatec-fn/xdasinh.f, cruft/slatec-fn/xdatanh.f, cruft/slatec-fn/xdbetai.f, cruft/slatec-fn/xderf.f, cruft/slatec-fn/xderfc.f, cruft/slatec-fn/xdgami.f, cruft/slatec-fn/xdgamit.f, cruft/slatec-fn/xdgamma.f, cruft/slatec-fn/xerf.f, cruft/slatec-fn/xerfc.f, cruft/slatec-fn/xgamma.f, cruft/slatec-fn/xgmainc.f, cruft/slatec-fn/xsgmainc.f: Moved from top-level libcruft to cruft directory below liboctave. * numeric/CmplxAEPBAL.cc, numeric/CmplxAEPBAL.h, numeric/CmplxCHOL.cc, numeric/CmplxCHOL.h, numeric/CmplxGEPBAL.cc, numeric/CmplxGEPBAL.h, numeric/CmplxHESS.cc, numeric/CmplxHESS.h, numeric/CmplxLU.cc, numeric/CmplxLU.h, numeric/CmplxQR.cc, numeric/CmplxQR.h, numeric/CmplxQRP.cc, numeric/CmplxQRP.h, numeric/CmplxSCHUR.cc, numeric/CmplxSCHUR.h, numeric/CmplxSVD.cc, numeric/CmplxSVD.h, numeric/CollocWt.cc, numeric/CollocWt.h, numeric/DAE.h, numeric/DAEFunc.h, numeric/DAERT.h, numeric/DAERTFunc.h, numeric/DASPK-opts.in, numeric/DASPK.cc, numeric/DASPK.h, numeric/DASRT-opts.in, numeric/DASRT.cc, numeric/DASRT.h, numeric/DASSL-opts.in, numeric/DASSL.cc, numeric/DASSL.h, numeric/DET.h, numeric/EIG.cc, numeric/EIG.h, numeric/LSODE-opts.in, numeric/LSODE.cc, numeric/LSODE.h, numeric/ODE.h, numeric/ODEFunc.h, numeric/ODES.cc, numeric/ODES.h, numeric/ODESFunc.h, numeric/Quad-opts.in, numeric/Quad.cc, numeric/Quad.h, numeric/SparseCmplxCHOL.cc, numeric/SparseCmplxCHOL.h, numeric/SparseCmplxLU.cc, numeric/SparseCmplxLU.h, numeric/SparseCmplxQR.cc, numeric/SparseCmplxQR.h, numeric/SparseQR.cc, numeric/SparseQR.h, numeric/SparsedbleCHOL.cc, numeric/SparsedbleCHOL.h, numeric/SparsedbleLU.cc, numeric/SparsedbleLU.h, numeric/base-aepbal.h, numeric/base-dae.h, numeric/base-de.h, numeric/base-lu.cc, numeric/base-lu.h, numeric/base-min.h, numeric/base-qr.cc, numeric/base-qr.h, numeric/bsxfun-decl.h, numeric/bsxfun-defs.cc, numeric/bsxfun.h, numeric/dbleAEPBAL.cc, numeric/dbleAEPBAL.h, numeric/dbleCHOL.cc, numeric/dbleCHOL.h, numeric/dbleGEPBAL.cc, numeric/dbleGEPBAL.h, numeric/dbleHESS.cc, numeric/dbleHESS.h, numeric/dbleLU.cc, numeric/dbleLU.h, numeric/dbleQR.cc, numeric/dbleQR.h, numeric/dbleQRP.cc, numeric/dbleQRP.h, numeric/dbleSCHUR.cc, numeric/dbleSCHUR.h, numeric/dbleSVD.cc, numeric/dbleSVD.h, numeric/eigs-base.cc, numeric/fCmplxAEPBAL.cc, numeric/fCmplxAEPBAL.h, numeric/fCmplxCHOL.cc, numeric/fCmplxCHOL.h, numeric/fCmplxGEPBAL.cc, numeric/fCmplxGEPBAL.h, numeric/fCmplxHESS.cc, numeric/fCmplxHESS.h, numeric/fCmplxLU.cc, numeric/fCmplxLU.h, numeric/fCmplxQR.cc, numeric/fCmplxQR.h, numeric/fCmplxQRP.cc, numeric/fCmplxQRP.h, numeric/fCmplxSCHUR.cc, numeric/fCmplxSCHUR.h, numeric/fCmplxSVD.cc, numeric/fCmplxSVD.h, numeric/fEIG.cc, numeric/fEIG.h, numeric/floatAEPBAL.cc, numeric/floatAEPBAL.h, numeric/floatCHOL.cc, numeric/floatCHOL.h, numeric/floatGEPBAL.cc, numeric/floatGEPBAL.h, numeric/floatHESS.cc, numeric/floatHESS.h, numeric/floatLU.cc, numeric/floatLU.h, numeric/floatQR.cc, numeric/floatQR.h, numeric/floatQRP.cc, numeric/floatQRP.h, numeric/floatSCHUR.cc, numeric/floatSCHUR.h, numeric/floatSVD.cc, numeric/floatSVD.h, numeric/lo-mappers.cc, numeric/lo-mappers.h, numeric/lo-specfun.cc, numeric/lo-specfun.h, numeric/module.mk, numeric/oct-convn.cc, numeric/oct-convn.h, numeric/oct-fftw.cc, numeric/oct-fftw.h, numeric/oct-norm.cc, numeric/oct-norm.h, numeric/oct-rand.cc, numeric/oct-rand.h, numeric/oct-spparms.cc, numeric/oct-spparms.h, numeric/randgamma.c, numeric/randgamma.h, numeric/randmtzig.c, numeric/randmtzig.h, numeric/randpoisson.c, numeric/randpoisson.h, numeric/sparse-base-chol.cc, numeric/sparse-base-chol.h, numeric/sparse-base-lu.cc, numeric/sparse-base-lu.h, numeric/sparse-dmsolve.cc: Moved from liboctave dir to numeric subdirectory. * operators/Sparse-diag-op-defs.h, operators/Sparse-op-defs.h, operators/Sparse-perm-op-defs.h, operators/config-ops.sh, operators/mk-ops.awk, operators/module.mk, operators/mx-base.h, operators/mx-defs.h, operators/mx-ext.h, operators/mx-inlines.cc, operators/mx-op-decl.h, operators/mx-op-defs.h, operators/mx-ops, operators/sparse-mk-ops.awk, operators/sparse-mx-ops, operators/vx-ops: Moved from liboctave dir to operators subdirectory. * system/dir-ops.cc, system/dir-ops.h, system/file-ops.cc, system/file-ops.h, system/file-stat.cc, system/file-stat.h, system/lo-sysdep.cc, system/lo-sysdep.h, system/mach-info.cc, system/mach-info.h, system/module.mk, system/oct-env.cc, system/oct-env.h, system/oct-group.cc, system/oct-group.h, system/oct-openmp.h, system/oct-passwd.cc, system/oct-passwd.h, system/oct-syscalls.cc, system/oct-syscalls.h, system/oct-time.cc, system/oct-time.h, system/oct-uname.cc, system/oct-uname.h, system/pathlen.h, system/sysdir.h, system/syswait.h, system/tempnam.c, system/tempname.c: Moved from liboctave dir to system subdirectory. * util/base-list.h, util/byte-swap.h, util/caseless-str.h, util/cmd-edit.cc, util/cmd-edit.h, util/cmd-hist.cc, util/cmd-hist.h, util/data-conv.cc, util/data-conv.h, util/f2c-main.c, util/functor.h, util/glob-match.cc, util/glob-match.h, util/kpse.cc, util/lo-array-gripes.cc, util/lo-array-gripes.h, util/lo-cieee.c, util/lo-cutils.c, util/lo-cutils.h, util/lo-ieee.cc, util/lo-ieee.h, util/lo-macros.h, util/lo-math.h, util/lo-traits.h, util/lo-utils.cc, util/lo-utils.h, util/module.mk, util/oct-alloc.cc, util/oct-alloc.h, util/oct-base64.cc, util/oct-base64.h, util/oct-binmap.h, util/oct-cmplx.h, util/oct-glob.cc, util/oct-glob.h, util/oct-inttypes.cc, util/oct-inttypes.h, util/oct-locbuf.cc, util/oct-locbuf.h, util/oct-md5.cc, util/oct-md5.h, util/oct-mem.h, util/oct-mutex.cc, util/oct-mutex.h, util/oct-refcount.h, util/oct-rl-edit.c, util/oct-rl-edit.h, util/oct-rl-hist.c, util/oct-rl-hist.h, util/oct-shlib.cc, util/oct-shlib.h, util/oct-sort.cc, util/oct-sort.h, util/oct-sparse.h, util/pathsearch.cc, util/pathsearch.h, util/regexp.cc, util/regexp.h, util/singleton-cleanup.cc, util/singleton-cleanup.h, util/sparse-sort.cc, util/sparse-sort.h, util/sparse-util.cc, util/sparse-util.h, util/statdefs.h, util/str-vec.cc, util/str-vec.h, util/sun-utils.h: Moved from liboctave dir to util subdirectory. * Makefile.am: Eliminate reference to top-level liboctave directory. * autogen.sh: cd to new liboctave/operators directory to run config-ops.sh. * build-aux/common.mk: Eliminate LIBCRUFT references. * configure.ac: Eliminate libcruft top-level references. Switch test programs to find files in liboctave/cruft subdirectory. * OctaveFAQ.texi, install.txi, mkoctfile.1: Eliminate references to libcruft in docs. * libgui/src/Makefile.am, libinterp/Makefile.am, src/Makefile.am: Update include file locations. Stop linking against libcruft. * libinterp/corefcn/module.mk: Update location of OPT_INC files which are now in numeric/ subdirectory. * libinterp/dldfcn/config-module.awk: Stop linking against libcruft. * libinterp/interpfcn/toplev.cc: Remove reference to LIBCRUFT. * libinterp/link-deps.mk, liboctave/link-deps.mk: Add GNULIB_LINK_DEPS to link dependencies. * libinterp/oct-conf.in.h: Remove reference to OCTAVE_CONF_LIBCRUFT. * liboctave/Makefile.am: Overhaul to use convenience libraries in subdirectories. * scripts/miscellaneous/mkoctfile.m: Eliminate reference to LIBCRUFT. * src/mkoctfile.in.cc, src/mkoctfile.in.sh: Stop linking againt libcruft. Eliminate references to LIBCRUFT.
author Rik <rik@octave.org>
date Fri, 31 Aug 2012 20:00:20 -0700
parents liboctave/Array.cc@9020dddc925a
children 7ed397c8ca68
comparison
equal deleted inserted replaced
15270:6615a46d90ec 15271:648dabbb4c6b
1 // Template array classes
2 /*
3
4 Copyright (C) 1993-2012 John W. Eaton
5 Copyright (C) 2008-2009 Jaroslav Hajek
6 Copyright (C) 2009 VZLU Prague
7
8 This file is part of Octave.
9
10 Octave is free software; you can redistribute it and/or modify it
11 under the terms of the GNU General Public License as published by the
12 Free Software Foundation; either version 3 of the License, or (at your
13 option) any later version.
14
15 Octave is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 for more details.
19
20 You should have received a copy of the GNU General Public License
21 along with Octave; see the file COPYING. If not, see
22 <http://www.gnu.org/licenses/>.
23
24 */
25
26 #ifdef HAVE_CONFIG_H
27 #include <config.h>
28 #endif
29
30 #include <cassert>
31
32 #include <iostream>
33 #include <sstream>
34 #include <vector>
35 #include <algorithm>
36 #include <new>
37
38 #include "Array.h"
39 #include "Array-util.h"
40 #include "idx-vector.h"
41 #include "lo-error.h"
42 #include "oct-locbuf.h"
43
44 // One dimensional array class. Handles the reference counting for
45 // all the derived classes.
46
47 template <class T>
48 Array<T>::Array (const Array<T>& a, const dim_vector& dv)
49 : dimensions (dv), rep (a.rep),
50 slice_data (a.slice_data), slice_len (a.slice_len)
51 {
52 if (dimensions.safe_numel () != a.numel ())
53 {
54 std::string dimensions_str = a.dimensions.str ();
55 std::string new_dims_str = dimensions.str ();
56
57 (*current_liboctave_error_handler)
58 ("reshape: can't reshape %s array to %s array",
59 dimensions_str.c_str (), new_dims_str.c_str ());
60 }
61
62 // This goes here because if an exception is thrown by the above,
63 // destructor will be never called.
64 rep->count++;
65 dimensions.chop_trailing_singletons ();
66 }
67
68 template <class T>
69 void
70 Array<T>::fill (const T& val)
71 {
72 if (rep->count > 1)
73 {
74 --rep->count;
75 rep = new ArrayRep (length (), val);
76 slice_data = rep->data;
77 }
78 else
79 fill_or_memset (slice_len, val, slice_data);
80 }
81
82 template <class T>
83 void
84 Array<T>::clear (void)
85 {
86 if (--rep->count == 0)
87 delete rep;
88
89 rep = nil_rep ();
90 rep->count++;
91 slice_data = rep->data;
92 slice_len = rep->len;
93
94 dimensions = dim_vector ();
95 }
96
97 template <class T>
98 void
99 Array<T>::clear (const dim_vector& dv)
100 {
101 if (--rep->count == 0)
102 delete rep;
103
104 rep = new ArrayRep (dv.safe_numel ());
105 slice_data = rep->data;
106 slice_len = rep->len;
107
108 dimensions = dv;
109 dimensions.chop_trailing_singletons ();
110 }
111
112 template <class T>
113 Array<T>
114 Array<T>::squeeze (void) const
115 {
116 Array<T> retval = *this;
117
118 if (ndims () > 2)
119 {
120 bool dims_changed = false;
121
122 dim_vector new_dimensions = dimensions;
123
124 int k = 0;
125
126 for (int i = 0; i < ndims (); i++)
127 {
128 if (dimensions(i) == 1)
129 dims_changed = true;
130 else
131 new_dimensions(k++) = dimensions(i);
132 }
133
134 if (dims_changed)
135 {
136 switch (k)
137 {
138 case 0:
139 new_dimensions = dim_vector (1, 1);
140 break;
141
142 case 1:
143 {
144 octave_idx_type tmp = new_dimensions(0);
145
146 new_dimensions.resize (2);
147
148 new_dimensions(0) = tmp;
149 new_dimensions(1) = 1;
150 }
151 break;
152
153 default:
154 new_dimensions.resize (k);
155 break;
156 }
157 }
158
159 retval = Array<T> (*this, new_dimensions);
160 }
161
162 return retval;
163 }
164
165 template <class T>
166 octave_idx_type
167 Array<T>::compute_index (octave_idx_type i, octave_idx_type j) const
168 {
169 return ::compute_index (i, j, dimensions);
170 }
171
172 template <class T>
173 octave_idx_type
174 Array<T>::compute_index (octave_idx_type i, octave_idx_type j, octave_idx_type k) const
175 {
176 return ::compute_index (i, j, k, dimensions);
177 }
178
179 template <class T>
180 octave_idx_type
181 Array<T>::compute_index (const Array<octave_idx_type>& ra_idx) const
182 {
183 return ::compute_index (ra_idx, dimensions);
184 }
185
186 template <class T>
187 T&
188 Array<T>::checkelem (octave_idx_type n)
189 {
190 // Do checks directly to avoid recomputing slice_len.
191 if (n < 0)
192 gripe_invalid_index ();
193 if (n >= slice_len)
194 gripe_index_out_of_range (1, 1, n+1, slice_len);
195
196 return elem (n);
197 }
198
199 template <class T>
200 T&
201 Array<T>::checkelem (octave_idx_type i, octave_idx_type j)
202 {
203 return elem (compute_index (i, j));
204 }
205
206 template <class T>
207 T&
208 Array<T>::checkelem (octave_idx_type i, octave_idx_type j, octave_idx_type k)
209 {
210 return elem (compute_index (i, j, k));
211 }
212
213 template <class T>
214 T&
215 Array<T>::checkelem (const Array<octave_idx_type>& ra_idx)
216 {
217 return elem (compute_index (ra_idx));
218 }
219
220 template <class T>
221 typename Array<T>::crefT
222 Array<T>::checkelem (octave_idx_type n) const
223 {
224 // Do checks directly to avoid recomputing slice_len.
225 if (n < 0)
226 gripe_invalid_index ();
227 if (n >= slice_len)
228 gripe_index_out_of_range (1, 1, n+1, slice_len);
229
230 return elem (n);
231 }
232
233 template <class T>
234 typename Array<T>::crefT
235 Array<T>::checkelem (octave_idx_type i, octave_idx_type j) const
236 {
237 return elem (compute_index (i, j));
238 }
239
240 template <class T>
241 typename Array<T>::crefT
242 Array<T>::checkelem (octave_idx_type i, octave_idx_type j, octave_idx_type k) const
243 {
244 return elem (compute_index (i, j, k));
245 }
246
247 template <class T>
248 typename Array<T>::crefT
249 Array<T>::checkelem (const Array<octave_idx_type>& ra_idx) const
250 {
251 return elem (compute_index (ra_idx));
252 }
253
254 template <class T>
255 Array<T>
256 Array<T>::column (octave_idx_type k) const
257 {
258 octave_idx_type r = dimensions(0);
259 #ifdef BOUNDS_CHECKING
260 if (k < 0 || k > dimensions.numel (1))
261 gripe_index_out_of_range (2, 2, k+1, dimensions.numel (1));
262 #endif
263
264 return Array<T> (*this, dim_vector (r, 1), k*r, k*r + r);
265 }
266
267 template <class T>
268 Array<T>
269 Array<T>::page (octave_idx_type k) const
270 {
271 octave_idx_type r = dimensions(0), c = dimensions (1), p = r*c;
272 #ifdef BOUNDS_CHECKING
273 if (k < 0 || k > dimensions.numel (2))
274 gripe_index_out_of_range (3, 3, k+1, dimensions.numel (2));
275 #endif
276
277 return Array<T> (*this, dim_vector (r, c), k*p, k*p + p);
278 }
279
280 template <class T>
281 Array<T>
282 Array<T>::linear_slice (octave_idx_type lo, octave_idx_type up) const
283 {
284 #ifdef BOUNDS_CHECKING
285 if (lo < 0)
286 gripe_index_out_of_range (1, 1, lo+1, numel ());
287 if (up > numel ())
288 gripe_index_out_of_range (1, 1, up, numel ());
289 #endif
290 if (up < lo) up = lo;
291 return Array<T> (*this, dim_vector (up - lo, 1), lo, up);
292 }
293
294 // Helper class for multi-d dimension permuting (generalized transpose).
295 class rec_permute_helper
296 {
297 // STRIDE occupies the last half of the space allocated for dim to
298 // avoid a double allocation.
299
300 int n;
301 int top;
302 octave_idx_type *dim;
303 octave_idx_type *stride;
304 bool use_blk;
305
306 public:
307 rec_permute_helper (const dim_vector& dv, const Array<octave_idx_type>& perm)
308
309 : n (dv.length ()), top (0), dim (new octave_idx_type [2*n]),
310 stride (dim + n), use_blk (false)
311 {
312 assert (n == perm.length ());
313
314 // Get cumulative dimensions.
315 OCTAVE_LOCAL_BUFFER (octave_idx_type, cdim, n+1);
316 cdim[0] = 1;
317 for (int i = 1; i < n+1; i++) cdim[i] = cdim[i-1] * dv(i-1);
318
319 // Setup the permuted strides.
320 for (int k = 0; k < n; k++)
321 {
322 int kk = perm(k);
323 dim[k] = dv(kk);
324 stride[k] = cdim[kk];
325 }
326
327 // Reduce contiguous runs.
328 for (int k = 1; k < n; k++)
329 {
330 if (stride[k] == stride[top]*dim[top])
331 dim[top] *= dim[k];
332 else
333 {
334 top++;
335 dim[top] = dim[k];
336 stride[top] = stride[k];
337 }
338 }
339
340 // Determine whether we can use block transposes.
341 use_blk = top >= 1 && stride[1] == 1 && stride[0] == dim[1];
342
343 }
344
345 ~rec_permute_helper (void) { delete [] dim; }
346
347 // Helper method for fast blocked transpose.
348 template <class T>
349 static T *
350 blk_trans (const T *src, T *dest, octave_idx_type nr, octave_idx_type nc)
351 {
352 static const octave_idx_type m = 8;
353 OCTAVE_LOCAL_BUFFER (T, blk, m*m);
354 for (octave_idx_type kr = 0; kr < nr; kr += m)
355 for (octave_idx_type kc = 0; kc < nc; kc += m)
356 {
357 octave_idx_type lr = std::min (m, nr - kr);
358 octave_idx_type lc = std::min (m, nc - kc);
359 if (lr == m && lc == m)
360 {
361 const T *ss = src + kc * nr + kr;
362 for (octave_idx_type j = 0; j < m; j++)
363 for (octave_idx_type i = 0; i < m; i++)
364 blk[j*m+i] = ss[j*nr + i];
365 T *dd = dest + kr * nc + kc;
366 for (octave_idx_type j = 0; j < m; j++)
367 for (octave_idx_type i = 0; i < m; i++)
368 dd[j*nc+i] = blk[i*m+j];
369 }
370 else
371 {
372 const T *ss = src + kc * nr + kr;
373 for (octave_idx_type j = 0; j < lc; j++)
374 for (octave_idx_type i = 0; i < lr; i++)
375 blk[j*m+i] = ss[j*nr + i];
376 T *dd = dest + kr * nc + kc;
377 for (octave_idx_type j = 0; j < lr; j++)
378 for (octave_idx_type i = 0; i < lc; i++)
379 dd[j*nc+i] = blk[i*m+j];
380 }
381 }
382
383 return dest + nr*nc;
384 }
385
386 private:
387
388 // Recursive N-d generalized transpose
389 template <class T>
390 T *do_permute (const T *src, T *dest, int lev) const
391 {
392 if (lev == 0)
393 {
394 octave_idx_type step = stride[0], len = dim[0];
395 if (step == 1)
396 {
397 copy_or_memcpy (len, src, dest);
398 dest += len;
399 }
400 else
401 {
402 for (octave_idx_type i = 0, j = 0; i < len; i++, j += step)
403 dest[i] = src[j];
404
405 dest += len;
406 }
407 }
408 else if (use_blk && lev == 1)
409 dest = blk_trans (src, dest, dim[1], dim[0]);
410 else
411 {
412 octave_idx_type step = stride[lev], len = dim[lev];
413 for (octave_idx_type i = 0, j = 0; i < len; i++, j+= step)
414 dest = do_permute (src + i * step, dest, lev-1);
415 }
416
417 return dest;
418 }
419
420 // No copying!
421
422 rec_permute_helper (const rec_permute_helper&);
423
424 rec_permute_helper& operator = (const rec_permute_helper&);
425
426 public:
427
428 template <class T>
429 void permute (const T *src, T *dest) const { do_permute (src, dest, top); }
430 };
431
432
433 template <class T>
434 Array<T>
435 Array<T>::permute (const Array<octave_idx_type>& perm_vec_arg, bool inv) const
436 {
437 Array<T> retval;
438
439 Array<octave_idx_type> perm_vec = perm_vec_arg;
440
441 dim_vector dv = dims ();
442
443 int perm_vec_len = perm_vec_arg.length ();
444
445 if (perm_vec_len < dv.length ())
446 (*current_liboctave_error_handler)
447 ("%s: invalid permutation vector", inv ? "ipermute" : "permute");
448
449 dim_vector dv_new = dim_vector::alloc (perm_vec_len);
450
451 // Append singleton dimensions as needed.
452 dv.resize (perm_vec_len, 1);
453
454 // Need this array to check for identical elements in permutation array.
455 OCTAVE_LOCAL_BUFFER_INIT (bool, checked, perm_vec_len, false);
456
457 bool identity = true;
458
459 // Find dimension vector of permuted array.
460 for (int i = 0; i < perm_vec_len; i++)
461 {
462 octave_idx_type perm_elt = perm_vec.elem (i);
463 if (perm_elt >= perm_vec_len || perm_elt < 0)
464 {
465 (*current_liboctave_error_handler)
466 ("%s: permutation vector contains an invalid element",
467 inv ? "ipermute" : "permute");
468
469 return retval;
470 }
471
472 if (checked[perm_elt])
473 {
474 (*current_liboctave_error_handler)
475 ("%s: permutation vector cannot contain identical elements",
476 inv ? "ipermute" : "permute");
477
478 return retval;
479 }
480 else
481 {
482 checked[perm_elt] = true;
483 identity = identity && perm_elt == i;
484 }
485 }
486
487 if (identity)
488 return *this;
489
490 if (inv)
491 {
492 for (int i = 0; i < perm_vec_len; i++)
493 perm_vec(perm_vec_arg(i)) = i;
494 }
495
496 for (int i = 0; i < perm_vec_len; i++)
497 dv_new(i) = dv(perm_vec(i));
498
499 retval = Array<T> (dv_new);
500
501 if (numel () > 0)
502 {
503 rec_permute_helper rh (dv, perm_vec);
504 rh.permute (data (), retval.fortran_vec ());
505 }
506
507 return retval;
508 }
509
510 // Helper class for multi-d index reduction and recursive
511 // indexing/indexed assignment. Rationale: we could avoid recursion
512 // using a state machine instead. However, using recursion is much
513 // more amenable to possible parallelization in the future.
514 // Also, the recursion solution is cleaner and more understandable.
515
516 class rec_index_helper
517 {
518 // CDIM occupies the last half of the space allocated for dim to
519 // avoid a double allocation.
520
521 int n;
522 int top;
523 octave_idx_type *dim;
524 octave_idx_type *cdim;
525 idx_vector *idx;
526
527 public:
528 rec_index_helper (const dim_vector& dv, const Array<idx_vector>& ia)
529 : n (ia.length ()), top (0), dim (new octave_idx_type [2*n]),
530 cdim (dim + n), idx (new idx_vector [n])
531 {
532 assert (n > 0 && (dv.length () == std::max (n, 2)));
533
534 dim[0] = dv(0);
535 cdim[0] = 1;
536 idx[0] = ia(0);
537
538 for (int i = 1; i < n; i++)
539 {
540 // Try reduction...
541 if (idx[top].maybe_reduce (dim[top], ia(i), dv(i)))
542 {
543 // Reduction successful, fold dimensions.
544 dim[top] *= dv(i);
545 }
546 else
547 {
548 // Unsuccessful, store index & cumulative dim.
549 top++;
550 idx[top] = ia(i);
551 dim[top] = dv(i);
552 cdim[top] = cdim[top-1] * dim[top-1];
553 }
554 }
555 }
556
557 ~rec_index_helper (void) { delete [] idx; delete [] dim; }
558
559 private:
560
561 // Recursive N-d indexing
562 template <class T>
563 T *do_index (const T *src, T *dest, int lev) const
564 {
565 if (lev == 0)
566 dest += idx[0].index (src, dim[0], dest);
567 else
568 {
569 octave_idx_type nn = idx[lev].length (dim[lev]), d = cdim[lev];
570 for (octave_idx_type i = 0; i < nn; i++)
571 dest = do_index (src + d*idx[lev].xelem (i), dest, lev-1);
572 }
573
574 return dest;
575 }
576
577 // Recursive N-d indexed assignment
578 template <class T>
579 const T *do_assign (const T *src, T *dest, int lev) const
580 {
581 if (lev == 0)
582 src += idx[0].assign (src, dim[0], dest);
583 else
584 {
585 octave_idx_type nn = idx[lev].length (dim[lev]), d = cdim[lev];
586 for (octave_idx_type i = 0; i < nn; i++)
587 src = do_assign (src, dest + d*idx[lev].xelem (i), lev-1);
588 }
589
590 return src;
591 }
592
593 // Recursive N-d indexed assignment
594 template <class T>
595 void do_fill (const T& val, T *dest, int lev) const
596 {
597 if (lev == 0)
598 idx[0].fill (val, dim[0], dest);
599 else
600 {
601 octave_idx_type nn = idx[lev].length (dim[lev]), d = cdim[lev];
602 for (octave_idx_type i = 0; i < nn; i++)
603 do_fill (val, dest + d*idx[lev].xelem (i), lev-1);
604 }
605 }
606
607 // No copying!
608
609 rec_index_helper (const rec_index_helper&);
610
611 rec_index_helper& operator = (const rec_index_helper&);
612
613 public:
614
615 template <class T>
616 void index (const T *src, T *dest) const { do_index (src, dest, top); }
617
618 template <class T>
619 void assign (const T *src, T *dest) const { do_assign (src, dest, top); }
620
621 template <class T>
622 void fill (const T& val, T *dest) const { do_fill (val, dest, top); }
623
624 bool is_cont_range (octave_idx_type& l,
625 octave_idx_type& u) const
626 {
627 return top == 0 && idx[0].is_cont_range (dim[0], l, u);
628 }
629 };
630
631 // Helper class for multi-d recursive resizing
632 // This handles resize () in an efficient manner, touching memory only
633 // once (apart from reinitialization)
634 class rec_resize_helper
635 {
636 octave_idx_type *cext;
637 octave_idx_type *sext;
638 octave_idx_type *dext;
639 int n;
640
641 public:
642 rec_resize_helper (const dim_vector& ndv, const dim_vector& odv)
643 : cext (0), sext (0), dext (0), n (0)
644 {
645 int l = ndv.length ();
646 assert (odv.length () == l);
647 octave_idx_type ld = 1;
648 int i = 0;
649 for (; i < l-1 && ndv(i) == odv(i); i++) ld *= ndv(i);
650 n = l - i;
651 cext = new octave_idx_type [3*n];
652 // Trick to avoid three allocations
653 sext = cext + n;
654 dext = sext + n;
655
656 octave_idx_type sld = ld, dld = ld;
657 for (int j = 0; j < n; j++)
658 {
659 cext[j] = std::min (ndv(i+j), odv(i+j));
660 sext[j] = sld *= odv(i+j);
661 dext[j] = dld *= ndv(i+j);
662 }
663 cext[0] *= ld;
664 }
665
666 ~rec_resize_helper (void) { delete [] cext; }
667
668 private:
669
670 // recursive resizing
671 template <class T>
672 void do_resize_fill (const T* src, T *dest, const T& rfv, int lev) const
673 {
674 if (lev == 0)
675 {
676 copy_or_memcpy (cext[0], src, dest);
677 fill_or_memset (dext[0] - cext[0], rfv, dest + cext[0]);
678 }
679 else
680 {
681 octave_idx_type sd = sext[lev-1], dd = dext[lev-1], k;
682 for (k = 0; k < cext[lev]; k++)
683 do_resize_fill (src + k * sd, dest + k * dd, rfv, lev - 1);
684
685 fill_or_memset (dext[lev] - k * dd, rfv, dest + k * dd);
686 }
687 }
688
689 // No copying!
690
691 rec_resize_helper (const rec_resize_helper&);
692
693 rec_resize_helper& operator = (const rec_resize_helper&);
694
695 public:
696
697 template <class T>
698 void resize_fill (const T* src, T *dest, const T& rfv) const
699 { do_resize_fill (src, dest, rfv, n-1); }
700 };
701
702 template <class T>
703 Array<T>
704 Array<T>::index (const idx_vector& i) const
705 {
706 octave_idx_type n = numel ();
707 Array<T> retval;
708
709 if (i.is_colon ())
710 {
711 // A(:) produces a shallow copy as a column vector.
712 retval = Array<T> (*this, dim_vector (n, 1));
713 }
714 else
715 {
716 if (i.extent (n) != n)
717 gripe_index_out_of_range (1, 1, i.extent (n), n); // throws
718
719 // FIXME -- this is the only place where orig_dimensions are used.
720 dim_vector rd = i.orig_dimensions ();
721 octave_idx_type il = i.length (n);
722
723 // FIXME -- this is for Matlab compatibility. Matlab 2007 given
724 //
725 // b = ones (3,1)
726 //
727 // yields the following:
728 //
729 // b(zeros (0,0)) gives []
730 // b(zeros (1,0)) gives zeros (0,1)
731 // b(zeros (0,1)) gives zeros (0,1)
732 // b(zeros (0,m)) gives zeros (0,m)
733 // b(zeros (m,0)) gives zeros (m,0)
734 // b(1:2) gives ones (2,1)
735 // b(ones (2)) gives ones (2) etc.
736 //
737 // As you can see, the behaviour is weird, but the tests end up pretty
738 // simple. Nah, I don't want to suggest that this is ad hoc :)
739
740 if (ndims () == 2 && n != 1 && rd.is_vector ())
741 {
742 if (columns () == 1)
743 rd = dim_vector (il, 1);
744 else if (rows () == 1)
745 rd = dim_vector (1, il);
746 }
747
748 octave_idx_type l, u;
749 if (il != 0 && i.is_cont_range (n, l, u))
750 // If suitable, produce a shallow slice.
751 retval = Array<T> (*this, rd, l, u);
752 else
753 {
754 // Don't use resize here to avoid useless initialization for POD
755 // types.
756 retval = Array<T> (rd);
757
758 if (il != 0)
759 i.index (data (), n, retval.fortran_vec ());
760 }
761 }
762
763 return retval;
764 }
765
766 template <class T>
767 Array<T>
768 Array<T>::index (const idx_vector& i, const idx_vector& j) const
769 {
770 // Get dimensions, allowing Fortran indexing in the 2nd dim.
771 dim_vector dv = dimensions.redim (2);
772 octave_idx_type r = dv(0), c = dv(1);
773 Array<T> retval;
774
775 if (i.is_colon () && j.is_colon ())
776 {
777 // A(:,:) produces a shallow copy.
778 retval = Array<T> (*this, dv);
779 }
780 else
781 {
782 if (i.extent (r) != r)
783 gripe_index_out_of_range (2, 1, i.extent (r), r); // throws
784 if (j.extent (c) != c)
785 gripe_index_out_of_range (2, 2, j.extent (c), c); // throws
786
787 octave_idx_type n = numel (), il = i.length (r), jl = j.length (c);
788
789 idx_vector ii (i);
790
791 if (ii.maybe_reduce (r, j, c))
792 {
793 octave_idx_type l, u;
794 if (ii.length () > 0 && ii.is_cont_range (n, l, u))
795 // If suitable, produce a shallow slice.
796 retval = Array<T> (*this, dim_vector (il, jl), l, u);
797 else
798 {
799 // Don't use resize here to avoid useless initialization for POD types.
800 retval = Array<T> (dim_vector (il, jl));
801
802 ii.index (data (), n, retval.fortran_vec ());
803 }
804 }
805 else
806 {
807 // Don't use resize here to avoid useless initialization for POD types.
808 retval = Array<T> (dim_vector (il, jl));
809
810 const T* src = data ();
811 T *dest = retval.fortran_vec ();
812
813 for (octave_idx_type k = 0; k < jl; k++)
814 dest += i.index (src + r * j.xelem (k), r, dest);
815 }
816 }
817
818 return retval;
819 }
820
821 template <class T>
822 Array<T>
823 Array<T>::index (const Array<idx_vector>& ia) const
824 {
825 int ial = ia.length ();
826 Array<T> retval;
827
828 // FIXME: is this dispatching necessary?
829 if (ial == 1)
830 retval = index (ia(0));
831 else if (ial == 2)
832 retval = index (ia(0), ia(1));
833 else if (ial > 0)
834 {
835 // Get dimensions, allowing Fortran indexing in the last dim.
836 dim_vector dv = dimensions.redim (ial);
837
838 // Check for out of bounds conditions.
839 bool all_colons = true;
840 for (int i = 0; i < ial; i++)
841 {
842 if (ia(i).extent (dv(i)) != dv(i))
843 gripe_index_out_of_range (ial, i+1, ia(i).extent (dv(i)), dv(i)); // throws
844
845 all_colons = all_colons && ia(i).is_colon ();
846 }
847
848
849 if (all_colons)
850 {
851 // A(:,:,...,:) produces a shallow copy.
852 dv.chop_trailing_singletons ();
853 retval = Array<T> (*this, dv);
854 }
855 else
856 {
857 // Form result dimensions.
858 dim_vector rdv = dim_vector::alloc (ial);
859 for (int i = 0; i < ial; i++) rdv(i) = ia(i).length (dv(i));
860 rdv.chop_trailing_singletons ();
861
862 // Prepare for recursive indexing
863 rec_index_helper rh (dv, ia);
864
865 octave_idx_type l, u;
866 if (rh.is_cont_range (l, u))
867 // If suitable, produce a shallow slice.
868 retval = Array<T> (*this, rdv, l, u);
869 else
870 {
871 // Don't use resize here to avoid useless initialization for POD types.
872 retval = Array<T> (rdv);
873
874 // Do it.
875 rh.index (data (), retval.fortran_vec ());
876 }
877 }
878 }
879
880 return retval;
881 }
882
883 // The default fill value. Override if you want a different one.
884
885 template <class T>
886 T
887 Array<T>::resize_fill_value (void) const
888 {
889 static T zero = T ();
890 return zero;
891 }
892
893 // Yes, we could do resize using index & assign. However, that would
894 // possibly involve a lot more memory traffic than we actually need.
895
896 template <class T>
897 void
898 Array<T>::resize1 (octave_idx_type n, const T& rfv)
899 {
900 if (n >= 0 && ndims () == 2)
901 {
902 dim_vector dv;
903 // This is driven by Matlab's behaviour of giving a *row* vector
904 // on some out-of-bounds assignments. Specifically, Matlab
905 // allows a(i) with out-of-bouds i when a is either of 0x0, 1x0,
906 // 1x1, 0xN, and gives a row vector in all cases (yes, even the
907 // last one, search me why). Giving a column vector would make
908 // much more sense (given the way trailing singleton dims are
909 // treated).
910 bool invalid = false;
911 if (rows () == 0 || rows () == 1)
912 dv = dim_vector (1, n);
913 else if (columns () == 1)
914 dv = dim_vector (n, 1);
915 else
916 invalid = true;
917
918 if (invalid)
919 gripe_invalid_resize ();
920 else
921 {
922 octave_idx_type nx = numel ();
923 if (n == nx - 1 && n > 0)
924 {
925 // Stack "pop" operation.
926 if (rep->count == 1)
927 slice_data[slice_len-1] = T ();
928 slice_len--;
929 dimensions = dv;
930 }
931 else if (n == nx + 1 && nx > 0)
932 {
933 // Stack "push" operation.
934 if (rep->count == 1 && slice_data + slice_len < rep->data + rep->len)
935 {
936 slice_data[slice_len++] = rfv;
937 dimensions = dv;
938 }
939 else
940 {
941 static const octave_idx_type max_stack_chunk = 1024;
942 octave_idx_type nn = n + std::min (nx, max_stack_chunk);
943 Array<T> tmp (Array<T> (dim_vector (nn, 1)), dv, 0, n);
944 T *dest = tmp.fortran_vec ();
945
946 copy_or_memcpy (nx, data (), dest);
947 dest[nx] = rfv;
948
949 *this = tmp;
950 }
951 }
952 else if (n != nx)
953 {
954 Array<T> tmp = Array<T> (dv);
955 T *dest = tmp.fortran_vec ();
956
957 octave_idx_type n0 = std::min (n, nx), n1 = n - n0;
958 copy_or_memcpy (n0, data (), dest);
959 fill_or_memset (n1, rfv, dest + n0);
960
961 *this = tmp;
962 }
963 }
964 }
965 else
966 gripe_invalid_resize ();
967 }
968
969 template <class T>
970 void
971 Array<T>::resize2 (octave_idx_type r, octave_idx_type c, const T& rfv)
972 {
973 if (r >= 0 && c >= 0 && ndims () == 2)
974 {
975 octave_idx_type rx = rows (), cx = columns ();
976 if (r != rx || c != cx)
977 {
978 Array<T> tmp = Array<T> (dim_vector (r, c));
979 T *dest = tmp.fortran_vec ();
980
981 octave_idx_type r0 = std::min (r, rx), r1 = r - r0;
982 octave_idx_type c0 = std::min (c, cx), c1 = c - c0;
983 const T *src = data ();
984 if (r == rx)
985 {
986 copy_or_memcpy (r * c0, src, dest);
987 dest += r * c0;
988 }
989 else
990 {
991 for (octave_idx_type k = 0; k < c0; k++)
992 {
993 copy_or_memcpy (r0, src, dest);
994 src += rx;
995 dest += r0;
996 fill_or_memset (r1, rfv, dest);
997 dest += r1;
998 }
999 }
1000
1001 fill_or_memset (r * c1, rfv, dest);
1002
1003 *this = tmp;
1004 }
1005 }
1006 else
1007 gripe_invalid_resize ();
1008
1009 }
1010
1011 template<class T>
1012 void
1013 Array<T>::resize (const dim_vector& dv, const T& rfv)
1014 {
1015 int dvl = dv.length ();
1016 if (dvl == 2)
1017 resize2 (dv(0), dv(1), rfv);
1018 else if (dimensions != dv)
1019 {
1020 if (dimensions.length () <= dvl && ! dv.any_neg ())
1021 {
1022 Array<T> tmp (dv);
1023 // Prepare for recursive resizing.
1024 rec_resize_helper rh (dv, dimensions.redim (dvl));
1025
1026 // Do it.
1027 rh.resize_fill (data (), tmp.fortran_vec (), rfv);
1028 *this = tmp;
1029 }
1030 else
1031 gripe_invalid_resize ();
1032 }
1033 }
1034
1035 template <class T>
1036 Array<T>
1037 Array<T>::index (const idx_vector& i, bool resize_ok, const T& rfv) const
1038 {
1039 Array<T> tmp = *this;
1040 if (resize_ok)
1041 {
1042 octave_idx_type n = numel (), nx = i.extent (n);
1043 if (n != nx)
1044 {
1045 if (i.is_scalar ())
1046 return Array<T> (dim_vector (1, 1), rfv);
1047 else
1048 tmp.resize1 (nx, rfv);
1049 }
1050
1051 if (tmp.numel () != nx)
1052 return Array<T> ();
1053 }
1054
1055 return tmp.index (i);
1056 }
1057
1058 template <class T>
1059 Array<T>
1060 Array<T>::index (const idx_vector& i, const idx_vector& j,
1061 bool resize_ok, const T& rfv) const
1062 {
1063 Array<T> tmp = *this;
1064 if (resize_ok)
1065 {
1066 dim_vector dv = dimensions.redim (2);
1067 octave_idx_type r = dv(0), c = dv(1);
1068 octave_idx_type rx = i.extent (r), cx = j.extent (c);
1069 if (r != rx || c != cx)
1070 {
1071 if (i.is_scalar () && j.is_scalar ())
1072 return Array<T> (dim_vector (1, 1), rfv);
1073 else
1074 tmp.resize2 (rx, cx, rfv);
1075 }
1076
1077 if (tmp.rows () != rx || tmp.columns () != cx)
1078 return Array<T> ();
1079 }
1080
1081 return tmp.index (i, j);
1082 }
1083
1084 template <class T>
1085 Array<T>
1086 Array<T>::index (const Array<idx_vector>& ia,
1087 bool resize_ok, const T& rfv) const
1088 {
1089 Array<T> tmp = *this;
1090 if (resize_ok)
1091 {
1092 int ial = ia.length ();
1093 dim_vector dv = dimensions.redim (ial);
1094 dim_vector dvx = dim_vector::alloc (ial);
1095 for (int i = 0; i < ial; i++) dvx(i) = ia(i).extent (dv (i));
1096 if (! (dvx == dv))
1097 {
1098 bool all_scalars = true;
1099 for (int i = 0; i < ial; i++)
1100 all_scalars = all_scalars && ia(i).is_scalar ();
1101 if (all_scalars)
1102 return Array<T> (dim_vector (1, 1), rfv);
1103 else
1104 tmp.resize (dvx, rfv);
1105 }
1106
1107 if (tmp.dimensions != dvx)
1108 return Array<T> ();
1109 }
1110
1111 return tmp.index (ia);
1112 }
1113
1114
1115 template <class T>
1116 void
1117 Array<T>::assign (const idx_vector& i, const Array<T>& rhs, const T& rfv)
1118 {
1119 octave_idx_type n = numel (), rhl = rhs.numel ();
1120
1121 if (rhl == 1 || i.length (n) == rhl)
1122 {
1123 octave_idx_type nx = i.extent (n);
1124 bool colon = i.is_colon_equiv (nx);
1125 // Try to resize first if necessary.
1126 if (nx != n)
1127 {
1128 // Optimize case A = []; A(1:n) = X with A empty.
1129 if (dimensions.zero_by_zero () && colon)
1130 {
1131 if (rhl == 1)
1132 *this = Array<T> (dim_vector (1, nx), rhs(0));
1133 else
1134 *this = Array<T> (rhs, dim_vector (1, nx));
1135 return;
1136 }
1137
1138 resize1 (nx, rfv);
1139 n = numel ();
1140 }
1141
1142 if (colon)
1143 {
1144 // A(:) = X makes a full fill or a shallow copy.
1145 if (rhl == 1)
1146 fill (rhs(0));
1147 else
1148 *this = rhs.reshape (dimensions);
1149 }
1150 else
1151 {
1152 if (rhl == 1)
1153 i.fill (rhs(0), n, fortran_vec ());
1154 else
1155 i.assign (rhs.data (), n, fortran_vec ());
1156 }
1157 }
1158 else
1159 gripe_invalid_assignment_size ();
1160 }
1161
1162 template <class T>
1163 void
1164 Array<T>::assign (const idx_vector& i, const idx_vector& j,
1165 const Array<T>& rhs, const T& rfv)
1166 {
1167 bool initial_dims_all_zero = dimensions.all_zero ();
1168
1169 // Get RHS extents, discarding singletons.
1170 dim_vector rhdv = rhs.dims ();
1171
1172 // Get LHS extents, allowing Fortran indexing in the second dim.
1173 dim_vector dv = dimensions.redim (2);
1174
1175 // Check for out-of-bounds and form resizing dimensions.
1176 dim_vector rdv;
1177
1178 // In the special when all dimensions are zero, colons are allowed
1179 // to inquire the shape of RHS. The rules are more obscure, so we
1180 // solve that elsewhere.
1181 if (initial_dims_all_zero)
1182 rdv = zero_dims_inquire (i, j, rhdv);
1183 else
1184 {
1185 rdv(0) = i.extent (dv(0));
1186 rdv(1) = j.extent (dv(1));
1187 }
1188
1189 bool isfill = rhs.numel () == 1;
1190 octave_idx_type il = i.length (rdv(0)), jl = j.length (rdv(1));
1191 rhdv.chop_all_singletons ();
1192 bool match = (isfill
1193 || (rhdv.length () == 2 && il == rhdv(0) && jl == rhdv(1)));
1194 match = match || (il == 1 && jl == rhdv(0) && rhdv(1) == 1);
1195
1196 if (match)
1197 {
1198 bool all_colons = (i.is_colon_equiv (rdv(0))
1199 && j.is_colon_equiv (rdv(1)));
1200 // Resize if requested.
1201 if (rdv != dv)
1202 {
1203 // Optimize case A = []; A(1:m, 1:n) = X
1204 if (dv.zero_by_zero () && all_colons)
1205 {
1206 if (isfill)
1207 *this = Array<T> (rdv, rhs(0));
1208 else
1209 *this = Array<T> (rhs, rdv);
1210 return;
1211 }
1212
1213 resize (rdv, rfv);
1214 dv = dimensions;
1215 }
1216
1217 if (all_colons)
1218 {
1219 // A(:,:) = X makes a full fill or a shallow copy
1220 if (isfill)
1221 fill (rhs(0));
1222 else
1223 *this = rhs.reshape (dimensions);
1224 }
1225 else
1226 {
1227 // The actual work.
1228 octave_idx_type n = numel (), r = dv (0), c = dv (1);
1229 idx_vector ii (i);
1230
1231 const T* src = rhs.data ();
1232 T *dest = fortran_vec ();
1233
1234 // Try reduction first.
1235 if (ii.maybe_reduce (r, j, c))
1236 {
1237 if (isfill)
1238 ii.fill (*src, n, dest);
1239 else
1240 ii.assign (src, n, dest);
1241 }
1242 else
1243 {
1244 if (isfill)
1245 {
1246 for (octave_idx_type k = 0; k < jl; k++)
1247 i.fill (*src, r, dest + r * j.xelem (k));
1248 }
1249 else
1250 {
1251 for (octave_idx_type k = 0; k < jl; k++)
1252 src += i.assign (src, r, dest + r * j.xelem (k));
1253 }
1254 }
1255 }
1256 }
1257 else
1258 gripe_assignment_dimension_mismatch ();
1259 }
1260
1261 template <class T>
1262 void
1263 Array<T>::assign (const Array<idx_vector>& ia,
1264 const Array<T>& rhs, const T& rfv)
1265 {
1266 int ial = ia.length ();
1267
1268 // FIXME: is this dispatching necessary / desirable?
1269 if (ial == 1)
1270 assign (ia(0), rhs, rfv);
1271 else if (ial == 2)
1272 assign (ia(0), ia(1), rhs, rfv);
1273 else if (ial > 0)
1274 {
1275 bool initial_dims_all_zero = dimensions.all_zero ();
1276
1277 // Get RHS extents, discarding singletons.
1278 dim_vector rhdv = rhs.dims ();
1279
1280 // Get LHS extents, allowing Fortran indexing in the second dim.
1281 dim_vector dv = dimensions.redim (ial);
1282
1283 // Get the extents forced by indexing.
1284 dim_vector rdv;
1285
1286 // In the special when all dimensions are zero, colons are
1287 // allowed to inquire the shape of RHS. The rules are more
1288 // obscure, so we solve that elsewhere.
1289 if (initial_dims_all_zero)
1290 rdv = zero_dims_inquire (ia, rhdv);
1291 else
1292 {
1293 rdv = dim_vector::alloc (ial);
1294 for (int i = 0; i < ial; i++)
1295 rdv(i) = ia(i).extent (dv(i));
1296 }
1297
1298 // Check whether LHS and RHS match, up to singleton dims.
1299 bool match = true, all_colons = true, isfill = rhs.numel () == 1;
1300
1301 rhdv.chop_all_singletons ();
1302 int j = 0, rhdvl = rhdv.length ();
1303 for (int i = 0; i < ial; i++)
1304 {
1305 all_colons = all_colons && ia(i).is_colon_equiv (rdv(i));
1306 octave_idx_type l = ia(i).length (rdv(i));
1307 if (l == 1) continue;
1308 match = match && j < rhdvl && l == rhdv(j++);
1309 }
1310
1311 match = match && (j == rhdvl || rhdv(j) == 1);
1312 match = match || isfill;
1313
1314 if (match)
1315 {
1316 // Resize first if necessary.
1317 if (rdv != dv)
1318 {
1319 // Optimize case A = []; A(1:m, 1:n) = X
1320 if (dv.zero_by_zero () && all_colons)
1321 {
1322 rdv.chop_trailing_singletons ();
1323 if (isfill)
1324 *this = Array<T> (rdv, rhs(0));
1325 else
1326 *this = Array<T> (rhs, rdv);
1327 return;
1328 }
1329
1330 resize (rdv, rfv);
1331 dv = rdv;
1332 }
1333
1334 if (all_colons)
1335 {
1336 // A(:,:,...,:) = X makes a full fill or a shallow copy.
1337 if (isfill)
1338 fill (rhs(0));
1339 else
1340 *this = rhs.reshape (dimensions);
1341 }
1342 else
1343 {
1344 // Do the actual work.
1345
1346 // Prepare for recursive indexing
1347 rec_index_helper rh (dv, ia);
1348
1349 // Do it.
1350 if (isfill)
1351 rh.fill (rhs(0), fortran_vec ());
1352 else
1353 rh.assign (rhs.data (), fortran_vec ());
1354 }
1355 }
1356 else
1357 gripe_assignment_dimension_mismatch ();
1358 }
1359 }
1360
1361 template <class T>
1362 void
1363 Array<T>::delete_elements (const idx_vector& i)
1364 {
1365 octave_idx_type n = numel ();
1366 if (i.is_colon ())
1367 {
1368 *this = Array<T> ();
1369 }
1370 else if (i.length (n) != 0)
1371 {
1372 if (i.extent (n) != n)
1373 gripe_del_index_out_of_range (true, i.extent (n), n);
1374
1375 octave_idx_type l, u;
1376 bool col_vec = ndims () == 2 && columns () == 1 && rows () != 1;
1377 if (i.is_scalar () && i(0) == n-1 && dimensions.is_vector ())
1378 {
1379 // Stack "pop" operation.
1380 resize1 (n-1);
1381 }
1382 else if (i.is_cont_range (n, l, u))
1383 {
1384 // Special case deleting a contiguous range.
1385 octave_idx_type m = n + l - u;
1386 Array<T> tmp (dim_vector (col_vec ? m : 1, !col_vec ? m : 1));
1387 const T *src = data ();
1388 T *dest = tmp.fortran_vec ();
1389 copy_or_memcpy (l, src, dest);
1390 copy_or_memcpy (n - u, src + u, dest + l);
1391 *this = tmp;
1392 }
1393 else
1394 {
1395 // Use index.
1396 *this = index (i.complement (n));
1397 }
1398 }
1399 }
1400
1401 template <class T>
1402 void
1403 Array<T>::delete_elements (int dim, const idx_vector& i)
1404 {
1405 if (dim < 0 || dim >= ndims ())
1406 {
1407 (*current_liboctave_error_handler)
1408 ("invalid dimension in delete_elements");
1409 return;
1410 }
1411
1412 octave_idx_type n = dimensions (dim);
1413 if (i.is_colon ())
1414 {
1415 *this = Array<T> ();
1416 }
1417 else if (i.length (n) != 0)
1418 {
1419 if (i.extent (n) != n)
1420 gripe_del_index_out_of_range (false, i.extent (n), n);
1421
1422 octave_idx_type l, u;
1423
1424 if (i.is_cont_range (n, l, u))
1425 {
1426 // Special case deleting a contiguous range.
1427 octave_idx_type nd = n + l - u, dl = 1, du = 1;
1428 dim_vector rdv = dimensions;
1429 rdv(dim) = nd;
1430 for (int k = 0; k < dim; k++) dl *= dimensions(k);
1431 for (int k = dim + 1; k < ndims (); k++) du *= dimensions(k);
1432
1433 // Special case deleting a contiguous range.
1434 Array<T> tmp = Array<T> (rdv);
1435 const T *src = data ();
1436 T *dest = tmp.fortran_vec ();
1437 l *= dl; u *= dl; n *= dl;
1438 for (octave_idx_type k = 0; k < du; k++)
1439 {
1440 copy_or_memcpy (l, src, dest);
1441 dest += l;
1442 copy_or_memcpy (n - u, src + u, dest);
1443 dest += n - u;
1444 src += n;
1445 }
1446
1447 *this = tmp;
1448 }
1449 else
1450 {
1451 // Use index.
1452 Array<idx_vector> ia (dim_vector (ndims (), 1), idx_vector::colon);
1453 ia (dim) = i.complement (n);
1454 *this = index (ia);
1455 }
1456 }
1457 }
1458
1459 template <class T>
1460 void
1461 Array<T>::delete_elements (const Array<idx_vector>& ia)
1462 {
1463 if (ia.length () == 1)
1464 delete_elements (ia(0));
1465 else
1466 {
1467 int len = ia.length (), k, dim = -1;
1468 for (k = 0; k < len; k++)
1469 {
1470 if (! ia(k).is_colon ())
1471 {
1472 if (dim < 0)
1473 dim = k;
1474 else
1475 break;
1476 }
1477 }
1478 if (dim < 0)
1479 {
1480 dim_vector dv = dimensions;
1481 dv(0) = 0;
1482 *this = Array<T> (dv);
1483 }
1484 else if (k == len)
1485 {
1486 delete_elements (dim, ia(dim));
1487 }
1488 else
1489 {
1490 (*current_liboctave_error_handler)
1491 ("a null assignment can only have one non-colon index");
1492 }
1493 }
1494
1495 }
1496
1497 template <class T>
1498 Array<T>&
1499 Array<T>::insert (const Array<T>& a, octave_idx_type r, octave_idx_type c)
1500 {
1501 idx_vector i (r, r + a.rows ());
1502 idx_vector j (c, c + a.columns ());
1503 if (ndims () == 2 && a.ndims () == 2)
1504 assign (i, j, a);
1505 else
1506 {
1507 Array<idx_vector> idx (dim_vector (a.ndims (), 1));
1508 idx(0) = i;
1509 idx(1) = j;
1510 for (int k = 2; k < a.ndims (); k++)
1511 idx(k) = idx_vector (0, a.dimensions(k));
1512 assign (idx, a);
1513 }
1514
1515 return *this;
1516 }
1517
1518 template <class T>
1519 Array<T>&
1520 Array<T>::insert (const Array<T>& a, const Array<octave_idx_type>& ra_idx)
1521 {
1522 octave_idx_type n = ra_idx.length ();
1523 Array<idx_vector> idx (dim_vector (n, 1));
1524 const dim_vector dva = a.dims ().redim (n);
1525 for (octave_idx_type k = 0; k < n; k++)
1526 idx(k) = idx_vector (ra_idx (k), ra_idx (k) + dva(k));
1527
1528 assign (idx, a);
1529
1530 return *this;
1531 }
1532
1533
1534 template <class T>
1535 Array<T>
1536 Array<T>::transpose (void) const
1537 {
1538 assert (ndims () == 2);
1539
1540 octave_idx_type nr = dim1 ();
1541 octave_idx_type nc = dim2 ();
1542
1543 if (nr >= 8 && nc >= 8)
1544 {
1545 Array<T> result (dim_vector (nc, nr));
1546
1547 // Reuse the implementation used for permuting.
1548
1549 rec_permute_helper::blk_trans (data (), result.fortran_vec (), nr, nc);
1550
1551 return result;
1552 }
1553 else if (nr > 1 && nc > 1)
1554 {
1555 Array<T> result (dim_vector (nc, nr));
1556
1557 for (octave_idx_type j = 0; j < nc; j++)
1558 for (octave_idx_type i = 0; i < nr; i++)
1559 result.xelem (j, i) = xelem (i, j);
1560
1561 return result;
1562 }
1563 else
1564 {
1565 // Fast transpose for vectors and empty matrices.
1566 return Array<T> (*this, dim_vector (nc, nr));
1567 }
1568 }
1569
1570 template <class T>
1571 static T
1572 no_op_fcn (const T& x)
1573 {
1574 return x;
1575 }
1576
1577 template <class T>
1578 Array<T>
1579 Array<T>::hermitian (T (*fcn) (const T&)) const
1580 {
1581 assert (ndims () == 2);
1582
1583 if (! fcn)
1584 fcn = no_op_fcn<T>;
1585
1586 octave_idx_type nr = dim1 ();
1587 octave_idx_type nc = dim2 ();
1588
1589 if (nr >= 8 && nc >= 8)
1590 {
1591 Array<T> result (dim_vector (nc, nr));
1592
1593 // Blocked transpose to attempt to avoid cache misses.
1594
1595 // Don't use OCTAVE_LOCAL_BUFFER here as it doesn't work with bool
1596 // on some compilers.
1597 T buf[64];
1598
1599 octave_idx_type ii = 0, jj;
1600 for (jj = 0; jj < (nc - 8 + 1); jj += 8)
1601 {
1602 for (ii = 0; ii < (nr - 8 + 1); ii += 8)
1603 {
1604 // Copy to buffer
1605 for (octave_idx_type j = jj, k = 0, idxj = jj * nr;
1606 j < jj + 8; j++, idxj += nr)
1607 for (octave_idx_type i = ii; i < ii + 8; i++)
1608 buf[k++] = xelem (i + idxj);
1609
1610 // Copy from buffer
1611 for (octave_idx_type i = ii, idxi = ii * nc; i < ii + 8;
1612 i++, idxi += nc)
1613 for (octave_idx_type j = jj, k = i - ii; j < jj + 8;
1614 j++, k+=8)
1615 result.xelem (j + idxi) = fcn (buf[k]);
1616 }
1617
1618 if (ii < nr)
1619 for (octave_idx_type j = jj; j < jj + 8; j++)
1620 for (octave_idx_type i = ii; i < nr; i++)
1621 result.xelem (j, i) = fcn (xelem (i, j));
1622 }
1623
1624 for (octave_idx_type j = jj; j < nc; j++)
1625 for (octave_idx_type i = 0; i < nr; i++)
1626 result.xelem (j, i) = fcn (xelem (i, j));
1627
1628 return result;
1629 }
1630 else
1631 {
1632 Array<T> result (dim_vector (nc, nr));
1633
1634 for (octave_idx_type j = 0; j < nc; j++)
1635 for (octave_idx_type i = 0; i < nr; i++)
1636 result.xelem (j, i) = fcn (xelem (i, j));
1637
1638 return result;
1639 }
1640 }
1641
1642 /*
1643
1644 %% Tranpose tests for matrices of the tile size and plus or minus a row
1645 %% and with four tiles.
1646
1647 %!shared m7, mt7, m8, mt8, m9, mt9
1648 %! m7 = reshape (1 : 7*8, 8, 7);
1649 %! mt7 = [1:8; 9:16; 17:24; 25:32; 33:40; 41:48; 49:56];
1650 %! m8 = reshape (1 : 8*8, 8, 8);
1651 %! mt8 = mt8 = [mt7; 57:64];
1652 %! m9 = reshape (1 : 9*8, 8, 9);
1653 %! mt9 = [mt8; 65:72];
1654
1655 %!assert (m7', mt7)
1656 %!assert ((1i*m7).', 1i * mt7)
1657 %!assert ((1i*m7)', conj (1i * mt7))
1658 %!assert (m8', mt8)
1659 %!assert ((1i*m8).', 1i * mt8)
1660 %!assert ((1i*m8)', conj (1i * mt8))
1661 %!assert (m9', mt9)
1662 %!assert ((1i*m9).', 1i * mt9)
1663 %!assert ((1i*m9)', conj (1i * mt9))
1664 %!assert ([m7, m8; m7, m8]', [mt7, mt7; mt8, mt8])
1665 %!assert ((1i*[m7, m8; m7, m8]).', 1i * [mt7, mt7; mt8, mt8])
1666 %!assert ((1i*[m7, m8; m7, m8])', conj (1i * [mt7, mt7; mt8, mt8]))
1667 %!assert ([m8, m8; m8, m8]', [mt8, mt8; mt8, mt8])
1668 %!assert ((1i*[m8, m8; m8, m8]).', 1i * [mt8, mt8; mt8, mt8])
1669 %!assert ((1i*[m8, m8; m8, m8])', conj (1i * [mt8, mt8; mt8, mt8]))
1670 %!assert ([m9, m8; m9, m8]', [mt9, mt9; mt8, mt8])
1671 %!assert ((1i*[m9, m8; m9, m8]).', 1i * [mt9, mt9; mt8, mt8])
1672 %!assert ((1i*[m9, m8; m9, m8])', conj (1i * [mt9, mt9; mt8, mt8]))
1673
1674 */
1675
1676 template <class T>
1677 T *
1678 Array<T>::fortran_vec (void)
1679 {
1680 make_unique ();
1681
1682 return slice_data;
1683 }
1684
1685 // Non-real types don't have NaNs.
1686 template <class T>
1687 inline bool
1688 sort_isnan (typename ref_param<T>::type)
1689 {
1690 return false;
1691 }
1692
1693 template <class T>
1694 Array<T>
1695 Array<T>::sort (int dim, sortmode mode) const
1696 {
1697 if (dim < 0)
1698 {
1699 (*current_liboctave_error_handler)
1700 ("sort: invalid dimension");
1701 return Array<T> ();
1702 }
1703
1704 Array<T> m (dims ());
1705
1706 dim_vector dv = m.dims ();
1707
1708 if (m.length () < 1)
1709 return m;
1710
1711 if (dim >= dv.length ())
1712 dv.resize (dim+1, 1);
1713
1714 octave_idx_type ns = dv(dim);
1715 octave_idx_type iter = dv.numel () / ns;
1716 octave_idx_type stride = 1;
1717
1718 for (int i = 0; i < dim; i++)
1719 stride *= dv(i);
1720
1721 T *v = m.fortran_vec ();
1722 const T *ov = data ();
1723
1724 octave_sort<T> lsort;
1725
1726 if (mode != UNSORTED)
1727 lsort.set_compare (mode);
1728 else
1729 return m;
1730
1731 if (stride == 1)
1732 {
1733 for (octave_idx_type j = 0; j < iter; j++)
1734 {
1735 // copy and partition out NaNs.
1736 // FIXME: impact on integer types noticeable?
1737 octave_idx_type kl = 0, ku = ns;
1738 for (octave_idx_type i = 0; i < ns; i++)
1739 {
1740 T tmp = ov[i];
1741 if (sort_isnan<T> (tmp))
1742 v[--ku] = tmp;
1743 else
1744 v[kl++] = tmp;
1745 }
1746
1747 // sort.
1748 lsort.sort (v, kl);
1749
1750 if (ku < ns)
1751 {
1752 // NaNs are in reverse order
1753 std::reverse (v + ku, v + ns);
1754 if (mode == DESCENDING)
1755 std::rotate (v, v + ku, v + ns);
1756 }
1757
1758 v += ns;
1759 ov += ns;
1760 }
1761 }
1762 else
1763 {
1764 OCTAVE_LOCAL_BUFFER (T, buf, ns);
1765
1766 for (octave_idx_type j = 0; j < iter; j++)
1767 {
1768 octave_idx_type offset = j;
1769 octave_idx_type offset2 = 0;
1770
1771 while (offset >= stride)
1772 {
1773 offset -= stride;
1774 offset2++;
1775 }
1776
1777 offset += offset2 * stride * ns;
1778
1779 // gather and partition out NaNs.
1780 // FIXME: impact on integer types noticeable?
1781 octave_idx_type kl = 0, ku = ns;
1782 for (octave_idx_type i = 0; i < ns; i++)
1783 {
1784 T tmp = ov[i*stride + offset];
1785 if (sort_isnan<T> (tmp))
1786 buf[--ku] = tmp;
1787 else
1788 buf[kl++] = tmp;
1789 }
1790
1791 // sort.
1792 lsort.sort (buf, kl);
1793
1794 if (ku < ns)
1795 {
1796 // NaNs are in reverse order
1797 std::reverse (buf + ku, buf + ns);
1798 if (mode == DESCENDING)
1799 std::rotate (buf, buf + ku, buf + ns);
1800 }
1801
1802 // scatter.
1803 for (octave_idx_type i = 0; i < ns; i++)
1804 v[i*stride + offset] = buf[i];
1805 }
1806 }
1807
1808 return m;
1809 }
1810
1811 template <class T>
1812 Array<T>
1813 Array<T>::sort (Array<octave_idx_type> &sidx, int dim,
1814 sortmode mode) const
1815 {
1816 if (dim < 0 || dim >= ndims ())
1817 {
1818 (*current_liboctave_error_handler)
1819 ("sort: invalid dimension");
1820 return Array<T> ();
1821 }
1822
1823 Array<T> m (dims ());
1824
1825 dim_vector dv = m.dims ();
1826
1827 if (m.length () < 1)
1828 {
1829 sidx = Array<octave_idx_type> (dv);
1830 return m;
1831 }
1832
1833 octave_idx_type ns = dv(dim);
1834 octave_idx_type iter = dv.numel () / ns;
1835 octave_idx_type stride = 1;
1836
1837 for (int i = 0; i < dim; i++)
1838 stride *= dv(i);
1839
1840 T *v = m.fortran_vec ();
1841 const T *ov = data ();
1842
1843 octave_sort<T> lsort;
1844
1845 sidx = Array<octave_idx_type> (dv);
1846 octave_idx_type *vi = sidx.fortran_vec ();
1847
1848 if (mode != UNSORTED)
1849 lsort.set_compare (mode);
1850 else
1851 return m;
1852
1853 if (stride == 1)
1854 {
1855 for (octave_idx_type j = 0; j < iter; j++)
1856 {
1857 // copy and partition out NaNs.
1858 // FIXME: impact on integer types noticeable?
1859 octave_idx_type kl = 0, ku = ns;
1860 for (octave_idx_type i = 0; i < ns; i++)
1861 {
1862 T tmp = ov[i];
1863 if (sort_isnan<T> (tmp))
1864 {
1865 --ku;
1866 v[ku] = tmp;
1867 vi[ku] = i;
1868 }
1869 else
1870 {
1871 v[kl] = tmp;
1872 vi[kl] = i;
1873 kl++;
1874 }
1875 }
1876
1877 // sort.
1878 lsort.sort (v, vi, kl);
1879
1880 if (ku < ns)
1881 {
1882 // NaNs are in reverse order
1883 std::reverse (v + ku, v + ns);
1884 std::reverse (vi + ku, vi + ns);
1885 if (mode == DESCENDING)
1886 {
1887 std::rotate (v, v + ku, v + ns);
1888 std::rotate (vi, vi + ku, vi + ns);
1889 }
1890 }
1891
1892 v += ns;
1893 vi += ns;
1894 ov += ns;
1895 }
1896 }
1897 else
1898 {
1899 OCTAVE_LOCAL_BUFFER (T, buf, ns);
1900 OCTAVE_LOCAL_BUFFER (octave_idx_type, bufi, ns);
1901
1902 for (octave_idx_type j = 0; j < iter; j++)
1903 {
1904 octave_idx_type offset = j;
1905 octave_idx_type offset2 = 0;
1906
1907 while (offset >= stride)
1908 {
1909 offset -= stride;
1910 offset2++;
1911 }
1912
1913 offset += offset2 * stride * ns;
1914
1915 // gather and partition out NaNs.
1916 // FIXME: impact on integer types noticeable?
1917 octave_idx_type kl = 0, ku = ns;
1918 for (octave_idx_type i = 0; i < ns; i++)
1919 {
1920 T tmp = ov[i*stride + offset];
1921 if (sort_isnan<T> (tmp))
1922 {
1923 --ku;
1924 buf[ku] = tmp;
1925 bufi[ku] = i;
1926 }
1927 else
1928 {
1929 buf[kl] = tmp;
1930 bufi[kl] = i;
1931 kl++;
1932 }
1933 }
1934
1935 // sort.
1936 lsort.sort (buf, bufi, kl);
1937
1938 if (ku < ns)
1939 {
1940 // NaNs are in reverse order
1941 std::reverse (buf + ku, buf + ns);
1942 std::reverse (bufi + ku, bufi + ns);
1943 if (mode == DESCENDING)
1944 {
1945 std::rotate (buf, buf + ku, buf + ns);
1946 std::rotate (bufi, bufi + ku, bufi + ns);
1947 }
1948 }
1949
1950 // scatter.
1951 for (octave_idx_type i = 0; i < ns; i++)
1952 v[i*stride + offset] = buf[i];
1953 for (octave_idx_type i = 0; i < ns; i++)
1954 vi[i*stride + offset] = bufi[i];
1955 }
1956 }
1957
1958 return m;
1959 }
1960
1961 template <class T>
1962 typename Array<T>::compare_fcn_type
1963 safe_comparator (sortmode mode, const Array<T>& /* a */,
1964 bool /* allow_chk */)
1965 {
1966 if (mode == ASCENDING)
1967 return octave_sort<T>::ascending_compare;
1968 else if (mode == DESCENDING)
1969 return octave_sort<T>::descending_compare;
1970 else
1971 return 0;
1972 }
1973
1974 template <class T>
1975 sortmode
1976 Array<T>::is_sorted (sortmode mode) const
1977 {
1978 octave_sort<T> lsort;
1979
1980 octave_idx_type n = numel ();
1981
1982 if (n <= 1)
1983 return mode ? mode : ASCENDING;
1984
1985 if (mode == UNSORTED)
1986 {
1987 // Auto-detect mode.
1988 compare_fcn_type compare
1989 = safe_comparator (ASCENDING, *this, false);
1990
1991 if (compare (elem (n-1), elem (0)))
1992 mode = DESCENDING;
1993 else
1994 mode = ASCENDING;
1995 }
1996
1997 if (mode != UNSORTED)
1998 {
1999 lsort.set_compare (safe_comparator (mode, *this, false));
2000
2001 if (! lsort.is_sorted (data (), n))
2002 mode = UNSORTED;
2003 }
2004
2005 return mode;
2006
2007 }
2008
2009 template <class T>
2010 Array<octave_idx_type>
2011 Array<T>::sort_rows_idx (sortmode mode) const
2012 {
2013 Array<octave_idx_type> idx;
2014
2015 octave_sort<T> lsort (safe_comparator (mode, *this, true));
2016
2017 octave_idx_type r = rows (), c = cols ();
2018
2019 idx = Array<octave_idx_type> (dim_vector (r, 1));
2020
2021 lsort.sort_rows (data (), idx.fortran_vec (), r, c);
2022
2023 return idx;
2024 }
2025
2026
2027 template <class T>
2028 sortmode
2029 Array<T>::is_sorted_rows (sortmode mode) const
2030 {
2031 octave_sort<T> lsort;
2032
2033 octave_idx_type r = rows (), c = cols ();
2034
2035 if (r <= 1 || c == 0)
2036 return mode ? mode : ASCENDING;
2037
2038 if (mode == UNSORTED)
2039 {
2040 // Auto-detect mode.
2041 compare_fcn_type compare
2042 = safe_comparator (ASCENDING, *this, false);
2043
2044 octave_idx_type i;
2045 for (i = 0; i < cols (); i++)
2046 {
2047 T l = elem (0, i), u = elem (rows () - 1, i);
2048 if (compare (l, u))
2049 {
2050 if (mode == DESCENDING)
2051 {
2052 mode = UNSORTED;
2053 break;
2054 }
2055 else
2056 mode = ASCENDING;
2057 }
2058 else if (compare (u, l))
2059 {
2060 if (mode == ASCENDING)
2061 {
2062 mode = UNSORTED;
2063 break;
2064 }
2065 else
2066 mode = DESCENDING;
2067 }
2068 }
2069 if (mode == UNSORTED && i == cols ())
2070 mode = ASCENDING;
2071 }
2072
2073 if (mode != UNSORTED)
2074 {
2075 lsort.set_compare (safe_comparator (mode, *this, false));
2076
2077 if (! lsort.is_sorted_rows (data (), r, c))
2078 mode = UNSORTED;
2079 }
2080
2081 return mode;
2082
2083 }
2084
2085 // Do a binary lookup in a sorted array.
2086 template <class T>
2087 octave_idx_type
2088 Array<T>::lookup (const T& value, sortmode mode) const
2089 {
2090 octave_idx_type n = numel ();
2091 octave_sort<T> lsort;
2092
2093 if (mode == UNSORTED)
2094 {
2095 // auto-detect mode
2096 if (n > 1 && lsort.descending_compare (elem (0), elem (n-1)))
2097 mode = DESCENDING;
2098 else
2099 mode = ASCENDING;
2100 }
2101
2102 lsort.set_compare (mode);
2103
2104 return lsort.lookup (data (), n, value);
2105 }
2106
2107 template <class T>
2108 Array<octave_idx_type>
2109 Array<T>::lookup (const Array<T>& values, sortmode mode) const
2110 {
2111 octave_idx_type n = numel (), nval = values.numel ();
2112 octave_sort<T> lsort;
2113 Array<octave_idx_type> idx (values.dims ());
2114
2115 if (mode == UNSORTED)
2116 {
2117 // auto-detect mode
2118 if (n > 1 && lsort.descending_compare (elem (0), elem (n-1)))
2119 mode = DESCENDING;
2120 else
2121 mode = ASCENDING;
2122 }
2123
2124 lsort.set_compare (mode);
2125
2126 // This determines the split ratio between the O(M*log2(N)) and O(M+N) algorithms.
2127 static const double ratio = 1.0;
2128 sortmode vmode = UNSORTED;
2129
2130 // Attempt the O(M+N) algorithm if M is large enough.
2131 if (nval > ratio * n / xlog2 (n + 1.0))
2132 {
2133 vmode = values.is_sorted ();
2134 // The table must not contain a NaN.
2135 if ((vmode == ASCENDING && sort_isnan<T> (values(nval-1)))
2136 || (vmode == DESCENDING && sort_isnan<T> (values(0))))
2137 vmode = UNSORTED;
2138 }
2139
2140 if (vmode != UNSORTED)
2141 lsort.lookup_sorted (data (), n, values.data (), nval,
2142 idx.fortran_vec (), vmode != mode);
2143 else
2144 lsort.lookup (data (), n, values.data (), nval, idx.fortran_vec ());
2145
2146 return idx;
2147 }
2148
2149 template <class T>
2150 octave_idx_type
2151 Array<T>::nnz (void) const
2152 {
2153 const T *src = data ();
2154 octave_idx_type nel = nelem (), retval = 0;
2155 const T zero = T ();
2156 for (octave_idx_type i = 0; i < nel; i++)
2157 if (src[i] != zero)
2158 retval++;
2159
2160 return retval;
2161 }
2162
2163 template <class T>
2164 Array<octave_idx_type>
2165 Array<T>::find (octave_idx_type n, bool backward) const
2166 {
2167 Array<octave_idx_type> retval;
2168 const T *src = data ();
2169 octave_idx_type nel = nelem ();
2170 const T zero = T ();
2171 if (n < 0 || n >= nel)
2172 {
2173 // We want all elements, which means we'll almost surely need
2174 // to resize. So count first, then allocate array of exact size.
2175 octave_idx_type cnt = 0;
2176 for (octave_idx_type i = 0; i < nel; i++)
2177 cnt += src[i] != zero;
2178
2179 retval.clear (cnt, 1);
2180 octave_idx_type *dest = retval.fortran_vec ();
2181 for (octave_idx_type i = 0; i < nel; i++)
2182 if (src[i] != zero) *dest++ = i;
2183 }
2184 else
2185 {
2186 // We want a fixed max number of elements, usually small. So be
2187 // optimistic, alloc the array in advance, and then resize if
2188 // needed.
2189 retval.clear (n, 1);
2190 if (backward)
2191 {
2192 // Do the search as a series of successive single-element searches.
2193 octave_idx_type k = 0, l = nel - 1;
2194 for (; k < n; k++)
2195 {
2196 for (;l >= 0 && src[l] == zero; l--) ;
2197 if (l >= 0)
2198 retval(k) = l--;
2199 else
2200 break;
2201 }
2202 if (k < n)
2203 retval.resize2 (k, 1);
2204 octave_idx_type *rdata = retval.fortran_vec ();
2205 std::reverse (rdata, rdata + k);
2206 }
2207 else
2208 {
2209 // Do the search as a series of successive single-element searches.
2210 octave_idx_type k = 0, l = 0;
2211 for (; k < n; k++)
2212 {
2213 for (;l != nel && src[l] == zero; l++) ;
2214 if (l != nel)
2215 retval(k) = l++;
2216 else
2217 break;
2218 }
2219 if (k < n)
2220 retval.resize2 (k, 1);
2221 }
2222 }
2223
2224 // Fixup return dimensions, for Matlab compatibility.
2225 // find (zeros (0,0)) -> zeros (0,0)
2226 // find (zeros (1,0)) -> zeros (1,0)
2227 // find (zeros (0,1)) -> zeros (0,1)
2228 // find (zeros (0,X)) -> zeros (0,1)
2229 // find (zeros (1,1)) -> zeros (0,0) !!!! WHY?
2230 // find (zeros (0,1,0)) -> zeros (0,0)
2231 // find (zeros (0,1,0,1)) -> zeros (0,0) etc
2232
2233 if ((numel () == 1 && retval.is_empty ())
2234 || (rows () == 0 && dims ().numel (1) == 0))
2235 retval.dimensions = dim_vector ();
2236 else if (rows () == 1 && ndims () == 2)
2237 retval.dimensions = dim_vector (1, retval.length ());
2238
2239 return retval;
2240 }
2241
2242 template <class T>
2243 Array<T>
2244 Array<T>::nth_element (const idx_vector& n, int dim) const
2245 {
2246 if (dim < 0)
2247 {
2248 (*current_liboctave_error_handler)
2249 ("nth_element: invalid dimension");
2250 return Array<T> ();
2251 }
2252
2253 dim_vector dv = dims ();
2254 if (dim >= dv.length ())
2255 dv.resize (dim+1, 1);
2256
2257 octave_idx_type ns = dv(dim);
2258
2259 octave_idx_type nn = n.length (ns);
2260
2261 dv(dim) = std::min (nn, ns);
2262 dv.chop_trailing_singletons ();
2263 dim = std::min (dv.length (), dim);
2264
2265 Array<T> m (dv);
2266
2267 if (m.numel () == 0)
2268 return m;
2269
2270 sortmode mode = UNSORTED;
2271 octave_idx_type lo = 0;
2272
2273 switch (n.idx_class ())
2274 {
2275 case idx_vector::class_scalar:
2276 mode = ASCENDING;
2277 lo = n(0);
2278 break;
2279 case idx_vector::class_range:
2280 {
2281 octave_idx_type inc = n.increment ();
2282 if (inc == 1)
2283 {
2284 mode = ASCENDING;
2285 lo = n(0);
2286 }
2287 else if (inc == -1)
2288 {
2289 mode = DESCENDING;
2290 lo = ns - 1 - n(0);
2291 }
2292 }
2293 default:
2294 break;
2295 }
2296
2297 if (mode == UNSORTED)
2298 {
2299 (*current_liboctave_error_handler)
2300 ("nth_element: n must be a scalar or a contiguous range");
2301 return Array<T> ();
2302 }
2303
2304 octave_idx_type up = lo + nn;
2305
2306 if (lo < 0 || up > ns)
2307 {
2308 (*current_liboctave_error_handler)
2309 ("nth_element: invalid element index");
2310 return Array<T> ();
2311 }
2312
2313 octave_idx_type iter = numel () / ns;
2314 octave_idx_type stride = 1;
2315
2316 for (int i = 0; i < dim; i++)
2317 stride *= dv(i);
2318
2319 T *v = m.fortran_vec ();
2320 const T *ov = data ();
2321
2322 OCTAVE_LOCAL_BUFFER (T, buf, ns);
2323
2324 octave_sort<T> lsort;
2325 lsort.set_compare (mode);
2326
2327 for (octave_idx_type j = 0; j < iter; j++)
2328 {
2329 octave_idx_type kl = 0, ku = ns;
2330
2331 if (stride == 1)
2332 {
2333 // copy without NaNs.
2334 // FIXME: impact on integer types noticeable?
2335 for (octave_idx_type i = 0; i < ns; i++)
2336 {
2337 T tmp = ov[i];
2338 if (sort_isnan<T> (tmp))
2339 buf[--ku] = tmp;
2340 else
2341 buf[kl++] = tmp;
2342 }
2343
2344 ov += ns;
2345 }
2346 else
2347 {
2348 octave_idx_type offset = j % stride;
2349 // copy without NaNs.
2350 // FIXME: impact on integer types noticeable?
2351 for (octave_idx_type i = 0; i < ns; i++)
2352 {
2353 T tmp = ov[offset + i*stride];
2354 if (sort_isnan<T> (tmp))
2355 buf[--ku] = tmp;
2356 else
2357 buf[kl++] = tmp;
2358 }
2359
2360 if (offset == stride-1)
2361 ov += ns*stride;
2362 }
2363
2364 if (ku == ns)
2365 lsort.nth_element (buf, ns, lo, up);
2366 else if (mode == ASCENDING)
2367 lsort.nth_element (buf, ku, lo, std::min (ku, up));
2368 else
2369 {
2370 octave_idx_type nnan = ns - ku;
2371 octave_idx_type zero = 0;
2372 lsort.nth_element (buf, ku, std::max (lo - nnan, zero),
2373 std::max (up - nnan, zero));
2374 std::rotate (buf, buf + ku, buf + ns);
2375 }
2376
2377 if (stride == 1)
2378 {
2379 for (octave_idx_type i = 0; i < nn; i++)
2380 v[i] = buf[lo + i];
2381
2382 v += nn;
2383 }
2384 else
2385 {
2386 octave_idx_type offset = j % stride;
2387 for (octave_idx_type i = 0; i < nn; i++)
2388 v[offset + stride * i] = buf[lo + i];
2389 if (offset == stride-1)
2390 v += nn*stride;
2391 }
2392 }
2393
2394 return m;
2395 }
2396
2397
2398 #define INSTANTIATE_ARRAY_SORT(T) template class OCTAVE_API octave_sort<T>;
2399
2400 #define NO_INSTANTIATE_ARRAY_SORT(T) \
2401 \
2402 template <> Array<T> \
2403 Array<T>::sort (int, sortmode) const { return *this; } \
2404 \
2405 template <> Array<T> \
2406 Array<T>::sort (Array<octave_idx_type> &sidx, int, sortmode) const \
2407 { sidx = Array<octave_idx_type> (); return *this; } \
2408 \
2409 template <> sortmode \
2410 Array<T>::is_sorted (sortmode) const \
2411 { return UNSORTED; } \
2412 \
2413 Array<T>::compare_fcn_type \
2414 safe_comparator (sortmode, const Array<T>&, bool) \
2415 { return 0; } \
2416 \
2417 template <> Array<octave_idx_type> \
2418 Array<T>::sort_rows_idx (sortmode) const \
2419 { return Array<octave_idx_type> (); } \
2420 \
2421 template <> sortmode \
2422 Array<T>::is_sorted_rows (sortmode) const \
2423 { return UNSORTED; } \
2424 \
2425 template <> octave_idx_type \
2426 Array<T>::lookup (T const &, sortmode) const \
2427 { return 0; } \
2428 template <> Array<octave_idx_type> \
2429 Array<T>::lookup (const Array<T>&, sortmode) const \
2430 { return Array<octave_idx_type> (); } \
2431 \
2432 template <> octave_idx_type \
2433 Array<T>::nnz (void) const\
2434 { return 0; } \
2435 template <> Array<octave_idx_type> \
2436 Array<T>::find (octave_idx_type, bool) const\
2437 { return Array<octave_idx_type> (); } \
2438 \
2439 template <> Array<T> \
2440 Array<T>::nth_element (const idx_vector&, int) const { return Array<T> (); } \
2441
2442
2443 template <class T>
2444 Array<T>
2445 Array<T>::diag (octave_idx_type k) const
2446 {
2447 dim_vector dv = dims ();
2448 octave_idx_type nd = dv.length ();
2449 Array<T> d;
2450
2451 if (nd > 2)
2452 (*current_liboctave_error_handler) ("Matrix must be 2-dimensional");
2453 else
2454 {
2455 octave_idx_type nnr = dv (0);
2456 octave_idx_type nnc = dv (1);
2457
2458 if (nnr == 0 && nnc == 0)
2459 ; // do nothing for empty matrix
2460 else if (nnr != 1 && nnc != 1)
2461 {
2462 // Extract diag from matrix
2463 if (k > 0)
2464 nnc -= k;
2465 else if (k < 0)
2466 nnr += k;
2467
2468 if (nnr > 0 && nnc > 0)
2469 {
2470 octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc;
2471
2472 d.resize (dim_vector (ndiag, 1));
2473
2474 if (k > 0)
2475 {
2476 for (octave_idx_type i = 0; i < ndiag; i++)
2477 d.xelem (i) = elem (i, i+k);
2478 }
2479 else if (k < 0)
2480 {
2481 for (octave_idx_type i = 0; i < ndiag; i++)
2482 d.xelem (i) = elem (i-k, i);
2483 }
2484 else
2485 {
2486 for (octave_idx_type i = 0; i < ndiag; i++)
2487 d.xelem (i) = elem (i, i);
2488 }
2489 }
2490 else
2491 (*current_liboctave_error_handler)
2492 ("diag: requested diagonal out of range");
2493 }
2494 else
2495 {
2496 // Create diag matrix from vector
2497 octave_idx_type roff = 0;
2498 octave_idx_type coff = 0;
2499 if (k > 0)
2500 {
2501 roff = 0;
2502 coff = k;
2503 }
2504 else if (k < 0)
2505 {
2506 roff = -k;
2507 coff = 0;
2508 }
2509
2510 if (nnr == 1)
2511 {
2512 octave_idx_type n = nnc + std::abs (k);
2513 d = Array<T> (dim_vector (n, n), resize_fill_value ());
2514
2515 for (octave_idx_type i = 0; i < nnc; i++)
2516 d.xelem (i+roff, i+coff) = elem (0, i);
2517 }
2518 else
2519 {
2520 octave_idx_type n = nnr + std::abs (k);
2521 d = Array<T> (dim_vector (n, n), resize_fill_value ());
2522
2523 for (octave_idx_type i = 0; i < nnr; i++)
2524 d.xelem (i+roff, i+coff) = elem (i, 0);
2525 }
2526 }
2527 }
2528
2529 return d;
2530 }
2531
2532 template <class T>
2533 Array<T>
2534 Array<T>::diag (octave_idx_type m, octave_idx_type n) const
2535 {
2536 Array<T> retval;
2537
2538 if (ndims () == 2 && (rows () == 1 || cols () == 1))
2539 {
2540 retval.resize (dim_vector (m, n), resize_fill_value ());
2541
2542 for (octave_idx_type i = 0; i < numel (); i++)
2543 retval.xelem (i, i) = xelem (i);
2544 }
2545 else
2546 (*current_liboctave_error_handler)
2547 ("cat: invalid dimension");
2548
2549 return retval;
2550 }
2551
2552 template <class T>
2553 Array<T>
2554 Array<T>::cat (int dim, octave_idx_type n, const Array<T> *array_list)
2555 {
2556 // Default concatenation.
2557 bool (dim_vector::*concat_rule) (const dim_vector&, int) = &dim_vector::concat;
2558
2559 if (dim == -1 || dim == -2)
2560 {
2561 concat_rule = &dim_vector::hvcat;
2562 dim = -dim - 1;
2563 }
2564 else if (dim < 0)
2565 (*current_liboctave_error_handler)
2566 ("cat: invalid dimension");
2567
2568 if (n == 1)
2569 return array_list[0];
2570 else if (n == 0)
2571 return Array<T> ();
2572
2573 // Special case:
2574 //
2575 // cat (dim, [], ..., [], A, ...)
2576 //
2577 // with dim > 2, A not 0x0, and at least three arguments to
2578 // concatenate is equivalent to
2579 //
2580 // cat (dim, A, ...)
2581 //
2582 // Note that this check must be performed here because for full-on
2583 // braindead Matlab compatibility, we need to have things like
2584 //
2585 // cat (3, [], [], A)
2586 //
2587 // succeed, but to have things like
2588 //
2589 // cat (3, cat (3, [], []), A)
2590 // cat (3, zeros (0, 0, 2), A)
2591 //
2592 // fail. See also bug report #31615.
2593
2594 octave_idx_type istart = 0;
2595
2596 if (n > 2 && dim > 1)
2597 {
2598 for (octave_idx_type i = 0; i < n; i++)
2599 {
2600 dim_vector dv = array_list[i].dims ();
2601
2602 if (dv.zero_by_zero ())
2603 istart++;
2604 else
2605 break;
2606 }
2607
2608 // Don't skip any initial aguments if they are all empty.
2609 if (istart >= n)
2610 istart = 0;
2611 }
2612
2613 dim_vector dv = array_list[istart++].dims ();
2614
2615 for (octave_idx_type i = istart; i < n; i++)
2616 if (! (dv.*concat_rule) (array_list[i].dims (), dim))
2617 (*current_liboctave_error_handler)
2618 ("cat: dimension mismatch");
2619
2620 Array<T> retval (dv);
2621
2622 if (retval.is_empty ())
2623 return retval;
2624
2625 int nidx = std::max (dv.length (), dim + 1);
2626 Array<idx_vector> idxa (dim_vector (nidx, 1), idx_vector::colon);
2627 octave_idx_type l = 0;
2628
2629 for (octave_idx_type i = 0; i < n; i++)
2630 {
2631 // NOTE: This takes some thinking, but no matter what the above rules
2632 // are, an empty array can always be skipped at this point, because
2633 // the result dimensions are already determined, and there is no way
2634 // an empty array may contribute a nonzero piece along the dimension
2635 // at this point, unless an empty array can be promoted to a non-empty
2636 // one (which makes no sense). I repeat, *no way*, think about it.
2637 if (array_list[i].is_empty ())
2638 continue;
2639
2640 octave_quit ();
2641
2642 octave_idx_type u;
2643 if (dim < array_list[i].ndims ())
2644 u = l + array_list[i].dims ()(dim);
2645 else
2646 u = l + 1;
2647
2648 idxa(dim) = idx_vector (l, u);
2649
2650 retval.assign (idxa, array_list[i]);
2651
2652 l = u;
2653 }
2654
2655 return retval;
2656 }
2657
2658 template <class T>
2659 void
2660 Array<T>::print_info (std::ostream& os, const std::string& prefix) const
2661 {
2662 os << prefix << "rep address: " << rep << '\n'
2663 << prefix << "rep->len: " << rep->len << '\n'
2664 << prefix << "rep->data: " << static_cast<void *> (rep->data) << '\n'
2665 << prefix << "rep->count: " << rep->count << '\n'
2666 << prefix << "slice_data: " << static_cast<void *> (slice_data) << '\n'
2667 << prefix << "slice_len: " << slice_len << '\n';
2668
2669 // 2D info:
2670 //
2671 // << pefix << "rows: " << rows () << "\n"
2672 // << prefix << "cols: " << cols () << "\n";
2673 }
2674
2675 template <class T>
2676 bool Array<T>::optimize_dimensions (const dim_vector& dv)
2677 {
2678 bool retval = dimensions == dv;
2679 if (retval)
2680 dimensions = dv;
2681
2682 return retval;
2683 }
2684
2685 template <class T>
2686 void Array<T>::instantiation_guard ()
2687 {
2688 // This guards against accidental implicit instantiations.
2689 // Array<T> instances should always be explicit and use INSTANTIATE_ARRAY.
2690 T::__xXxXx__ ();
2691 }
2692
2693 #define INSTANTIATE_ARRAY(T, API) \
2694 template <> void Array<T>::instantiation_guard () { } \
2695 template class API Array<T>
2696
2697 // FIXME: is this used?
2698
2699 template <class T>
2700 std::ostream&
2701 operator << (std::ostream& os, const Array<T>& a)
2702 {
2703 dim_vector a_dims = a.dims ();
2704
2705 int n_dims = a_dims.length ();
2706
2707 os << n_dims << "-dimensional array";
2708
2709 if (n_dims)
2710 os << " (" << a_dims.str () << ")";
2711
2712 os <<"\n\n";
2713
2714 if (n_dims)
2715 {
2716 os << "data:";
2717
2718 Array<octave_idx_type> ra_idx (dim_vector (n_dims, 1), 0);
2719
2720 // Number of times the first 2d-array is to be displayed.
2721
2722 octave_idx_type m = 1;
2723 for (int i = 2; i < n_dims; i++)
2724 m *= a_dims(i);
2725
2726 if (m == 1)
2727 {
2728 octave_idx_type rows = 0;
2729 octave_idx_type cols = 0;
2730
2731 switch (n_dims)
2732 {
2733 case 2:
2734 rows = a_dims(0);
2735 cols = a_dims(1);
2736
2737 for (octave_idx_type j = 0; j < rows; j++)
2738 {
2739 ra_idx(0) = j;
2740 for (octave_idx_type k = 0; k < cols; k++)
2741 {
2742 ra_idx(1) = k;
2743 os << " " << a.elem (ra_idx);
2744 }
2745 os << "\n";
2746 }
2747 break;
2748
2749 default:
2750 rows = a_dims(0);
2751
2752 for (octave_idx_type k = 0; k < rows; k++)
2753 {
2754 ra_idx(0) = k;
2755 os << " " << a.elem (ra_idx);
2756 }
2757 break;
2758 }
2759
2760 os << "\n";
2761 }
2762 else
2763 {
2764 octave_idx_type rows = a_dims(0);
2765 octave_idx_type cols = a_dims(1);
2766
2767 for (int i = 0; i < m; i++)
2768 {
2769 os << "\n(:,:,";
2770
2771 for (int j = 2; j < n_dims - 1; j++)
2772 os << ra_idx(j) + 1 << ",";
2773
2774 os << ra_idx(n_dims - 1) + 1 << ") = \n";
2775
2776 for (octave_idx_type j = 0; j < rows; j++)
2777 {
2778 ra_idx(0) = j;
2779
2780 for (octave_idx_type k = 0; k < cols; k++)
2781 {
2782 ra_idx(1) = k;
2783 os << " " << a.elem (ra_idx);
2784 }
2785
2786 os << "\n";
2787 }
2788
2789 os << "\n";
2790
2791 if (i != m - 1)
2792 increment_index (ra_idx, a_dims, 2);
2793 }
2794 }
2795 }
2796
2797 return os;
2798 }