Mercurial > octave
view liboctave/operators/Sparse-op-defs.h @ 21244:1473547f50f5
include octave-config.h in public header files
* mk-opts.pl, mkbuiltins, mk-ops.awk, sparse-mk-ops.awk:
Emit "#include octave-config.h" statement for generated header files.
* build-env.h, builtins.h, Cell.h, base-text-renderer.h,
c-file-ptr-stream.h, cdisplay.h, comment-list.h, data.h, debug.h,
defaults.in.h, defun-dld.h, defun-int.h, defun.h, dirfns.h, display.h,
dynamic-ld.h, error.h, errwarn.h, event-queue.h, file-io.h,
ft-text-renderer.h, gl-render.h, gl2ps-print.h, graphics.in.h,
gripes.h, help.h, hook-fcn.h, input.h, jit-ir.h, jit-typeinfo.h,
jit-util.h, load-path.h, load-save.h, ls-ascii-helper.h, ls-hdf5.h,
ls-mat-ascii.h, ls-mat4.h, ls-mat5.h, ls-oct-binary.h, ls-oct-text.h,
ls-utils.h, mex.h, mexproto.h, mxarray.in.h, oct-errno.h, oct-fstrm.h,
oct-handle.h, oct-hdf5-types.h, oct-hdf5.h, oct-hist.h, oct-iostrm.h,
oct-lvalue.h, oct-map.h, oct-obj.h, oct-opengl.h, oct-prcstrm.h,
oct-procbuf.h, oct-stdstrm.h, oct-stream.h, oct-strstrm.h, oct.h,
octave-default-image.h, octave-link.h, octave-preserve-stream-state.h,
pager.h, pr-output.h, procstream.h, profiler.h, pt-jit.h,
sighandlers.h, siglist.h, sparse-xdiv.h, sparse-xpow.h, symtab.h,
sysdep.h, text-renderer.h, toplev.h, txt-eng.h, utils.h, variables.h,
workspace-element.h, xdiv.h, xnorm.h, xpow.h, zfstream.h, oct-qhull.h,
ov-base-diag.h, ov-base-int.h, ov-base-mat.h, ov-base-scalar.h,
ov-base-sparse.h, ov-base.h, ov-bool-mat.h, ov-bool-sparse.h,
ov-bool.h, ov-builtin.h, ov-cell.h, ov-ch-mat.h, ov-class.h,
ov-classdef.h, ov-colon.h, ov-complex.h, ov-cs-list.h, ov-cx-diag.h,
ov-cx-mat.h, ov-cx-sparse.h, ov-dld-fcn.h, ov-fcn-handle.h,
ov-fcn-inline.h, ov-fcn.h, ov-float.h, ov-flt-complex.h,
ov-flt-cx-diag.h, ov-flt-cx-mat.h, ov-flt-re-diag.h, ov-flt-re-mat.h,
ov-int-traits.h, ov-int16.h, ov-int32.h, ov-int64.h, ov-int8.h,
ov-intx.h, ov-java.h, ov-lazy-idx.h, ov-mex-fcn.h, ov-null-mat.h,
ov-oncleanup.h, ov-perm.h, ov-range.h, ov-re-diag.h, ov-re-mat.h,
ov-re-sparse.h, ov-scalar.h, ov-str-mat.h, ov-struct.h,
ov-type-conv.h, ov-typeinfo.h, ov-uint16.h, ov-uint32.h, ov-uint64.h,
ov-uint8.h, ov-usr-fcn.h, ov.h, ovl.h, octave.h, op-int.h, ops.h,
options-usage.h, lex.h, parse.h, pt-all.h, pt-arg-list.h,
pt-array-list.h, pt-assign.h, pt-binop.h, pt-bp.h, pt-cbinop.h,
pt-cell.h, pt-check.h, pt-classdef.h, pt-cmd.h, pt-colon.h,
pt-const.h, pt-decl.h, pt-eval.h, pt-except.h, pt-exp.h,
pt-fcn-handle.h, pt-funcall.h, pt-id.h, pt-idx.h, pt-jump.h,
pt-loop.h, pt-mat.h, pt-misc.h, pt-pr-code.h, pt-select.h, pt-stmt.h,
pt-unop.h, pt-walk.h, pt.h, token.h, version.in.h, Array-util.h,
Array.h, CColVector.h, CDiagMatrix.h, CMatrix.h, CNDArray.h,
CRowVector.h, CSparse.h, DiagArray2.h, MArray.h, MDiagArray2.h,
MSparse.h, Matrix.h, MatrixType.h, PermMatrix.h, Range.h, Sparse.h,
boolMatrix.h, boolNDArray.h, boolSparse.h, chMatrix.h, chNDArray.h,
dColVector.h, dDiagMatrix.h, dMatrix.h, dNDArray.h, dRowVector.h,
dSparse.h, dim-vector.h, fCColVector.h, fCDiagMatrix.h, fCMatrix.h,
fCNDArray.h, fCRowVector.h, fColVector.h, fDiagMatrix.h, fMatrix.h,
fNDArray.h, fRowVector.h, idx-vector.h, int16NDArray.h,
int32NDArray.h, int64NDArray.h, int8NDArray.h, intNDArray.h,
uint16NDArray.h, uint32NDArray.h, uint64NDArray.h, uint8NDArray.h,
f77-fcn.h, lo-error.h, quit.h, CmplxAEPBAL.h, CmplxCHOL.h,
CmplxGEPBAL.h, CmplxHESS.h, CmplxLU.h, CmplxQR.h, CmplxQRP.h,
CmplxSCHUR.h, CmplxSVD.h, CollocWt.h, DAE.h, DAEFunc.h, DAERT.h,
DAERTFunc.h, DASPK.h, DASRT.h, DASSL.h, DET.h, EIG.h, LSODE.h, ODE.h,
ODEFunc.h, ODES.h, ODESFunc.h, Quad.h, base-aepbal.h, base-dae.h,
base-de.h, base-lu.h, base-min.h, base-qr.h, bsxfun-decl.h, bsxfun.h,
dbleAEPBAL.h, dbleCHOL.h, dbleGEPBAL.h, dbleHESS.h, dbleLU.h,
dbleQR.h, dbleQRP.h, dbleSCHUR.h, dbleSVD.h, eigs-base.h,
fCmplxAEPBAL.h, fCmplxCHOL.h, fCmplxGEPBAL.h, fCmplxHESS.h,
fCmplxLU.h, fCmplxQR.h, fCmplxQRP.h, fCmplxSCHUR.h, fCmplxSVD.h,
fEIG.h, floatAEPBAL.h, floatCHOL.h, floatGEPBAL.h, floatHESS.h,
floatLU.h, floatQR.h, floatQRP.h, floatSCHUR.h, floatSVD.h,
lo-mappers.h, lo-specfun.h, oct-convn.h, oct-fftw.h, oct-norm.h,
oct-rand.h, oct-spparms.h, randgamma.h, randmtzig.h, randpoisson.h,
sparse-chol.h, sparse-dmsolve.h, sparse-lu.h, sparse-qr.h,
Sparse-diag-op-defs.h, Sparse-op-decls.h, Sparse-op-defs.h,
Sparse-perm-op-defs.h, mx-base.h, mx-defs.h, mx-ext.h, mx-op-decl.h,
mx-op-defs.h, dir-ops.h, file-ops.h, file-stat.h, lo-sysdep.h,
mach-info.h, oct-env.h, oct-group.h, oct-openmp.h, oct-passwd.h,
oct-syscalls.h, oct-time.h, oct-uname.h, pathlen.h, sysdir.h,
syswait.h, action-container.h, base-list.h, byte-swap.h,
caseless-str.h, cmd-edit.h, cmd-hist.h, data-conv.h, functor.h,
glob-match.h, lo-array-errwarn.h, lo-array-gripes.h, lo-cutils.h,
lo-ieee.h, lo-macros.h, lo-math.h, lo-regexp.h, lo-traits.h,
lo-utils.h, oct-alloc.h, oct-base64.h, oct-binmap.h, oct-cmplx.h,
oct-glob.h, oct-inttypes.h, oct-locbuf.h, oct-mutex.h, oct-refcount.h,
oct-rl-edit.h, oct-rl-hist.h, oct-shlib.h, oct-sort.h, oct-sparse.h,
pathsearch.h, singleton-cleanup.h, sparse-sort.h, sparse-util.h,
statdefs.h, str-vec.h, sun-utils.h, unwind-prot.h, url-transfer.h:
Include octave-config.h.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Wed, 10 Feb 2016 14:25:53 -0500 |
| parents | 228b65504557 |
| children | 6c2fd62db1f7 |
line wrap: on
line source
/* Copyright (C) 2004-2015 David Bateman Copyright (C) 1998-2004 Andy Adler Copyright (C) 2008 Jaroslav Hajek This file is part of Octave. Octave is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. Octave is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with Octave; see the file COPYING. If not, see <http://www.gnu.org/licenses/>. */ #if ! defined (octave_Sparse_op_defs_h) #define octave_Sparse_op_defs_h 1 #include "octave-config.h" #include "Array-util.h" #include "oct-locbuf.h" #include "mx-inlines.cc" // sparse matrix by scalar operations. #define SPARSE_SMS_BIN_OP_1(R, F, OP, M, S) \ R \ F (const M& m, const S& s) \ { \ octave_idx_type nr = m.rows (); \ octave_idx_type nc = m.cols (); \ \ R r (nr, nc, (0.0 OP s)); \ \ for (octave_idx_type j = 0; j < nc; j++) \ for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ r.xelem (m.ridx (i), j) = m.data (i) OP s; \ return r; \ } #define SPARSE_SMS_BIN_OP_2(R, F, OP, M, S) \ R \ F (const M& m, const S& s) \ { \ octave_idx_type nr = m.rows (); \ octave_idx_type nc = m.cols (); \ octave_idx_type nz = m.nnz (); \ \ R r (nr, nc, nz); \ \ for (octave_idx_type i = 0; i < nz; i++) \ { \ r.xdata (i) = m.data (i) OP s; \ r.xridx (i) = m.ridx (i); \ } \ for (octave_idx_type i = 0; i < nc + 1; i++) \ r.xcidx (i) = m.cidx (i); \ \ r.maybe_compress (true); \ return r; \ } #define SPARSE_SMS_BIN_OPS(R1, R2, M, S) \ SPARSE_SMS_BIN_OP_1 (R1, operator +, +, M, S) \ SPARSE_SMS_BIN_OP_1 (R1, operator -, -, M, S) \ SPARSE_SMS_BIN_OP_2 (R2, operator *, *, M, S) \ SPARSE_SMS_BIN_OP_2 (R2, operator /, /, M, S) #define SPARSE_SMS_CMP_OP(F, OP, M, MZ, MC, S, SZ, SC) \ SparseBoolMatrix \ F (const M& m, const S& s) \ { \ octave_idx_type nr = m.rows (); \ octave_idx_type nc = m.cols (); \ SparseBoolMatrix r; \ \ if (MC (MZ) OP SC (s)) \ { \ r = SparseBoolMatrix (nr, nc, true); \ for (octave_idx_type j = 0; j < nc; j++) \ for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ if (! (MC (m.data (i)) OP SC (s))) \ r.data (m.ridx (i) + j * nr) = false; \ r.maybe_compress (true); \ } \ else \ { \ r = SparseBoolMatrix (nr, nc, m.nnz ()); \ r.cidx (0) = static_cast<octave_idx_type> (0); \ octave_idx_type nel = 0; \ for (octave_idx_type j = 0; j < nc; j++) \ { \ for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ if (MC (m.data (i)) OP SC (s)) \ { \ r.ridx (nel) = m.ridx (i); \ r.data (nel++) = true; \ } \ r.cidx (j + 1) = nel; \ } \ r.maybe_compress (false); \ } \ return r; \ } #define SPARSE_SMS_CMP_OPS(M, MZ, CM, S, SZ, CS) \ SPARSE_SMS_CMP_OP (mx_el_lt, <, M, MZ, , S, SZ, ) \ SPARSE_SMS_CMP_OP (mx_el_le, <=, M, MZ, , S, SZ, ) \ SPARSE_SMS_CMP_OP (mx_el_ge, >=, M, MZ, , S, SZ, ) \ SPARSE_SMS_CMP_OP (mx_el_gt, >, M, MZ, , S, SZ, ) \ SPARSE_SMS_CMP_OP (mx_el_eq, ==, M, MZ, , S, SZ, ) \ SPARSE_SMS_CMP_OP (mx_el_ne, !=, M, MZ, , S, SZ, ) #define SPARSE_SMS_EQNE_OPS(M, MZ, CM, S, SZ, CS) \ SPARSE_SMS_CMP_OP (mx_el_eq, ==, M, MZ, , S, SZ, ) \ SPARSE_SMS_CMP_OP (mx_el_ne, !=, M, MZ, , S, SZ, ) #define SPARSE_SMS_BOOL_OP(F, OP, M, S, LHS_ZERO, RHS_ZERO) \ SparseBoolMatrix \ F (const M& m, const S& s) \ { \ octave_idx_type nr = m.rows (); \ octave_idx_type nc = m.cols (); \ SparseBoolMatrix r; \ \ if (nr > 0 && nc > 0) \ { \ if (LHS_ZERO OP (s != RHS_ZERO)) \ { \ r = SparseBoolMatrix (nr, nc, true); \ for (octave_idx_type j = 0; j < nc; j++) \ for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ if (! ((m.data (i) != LHS_ZERO) OP (s != RHS_ZERO))) \ r.data (m.ridx (i) + j * nr) = false; \ r.maybe_compress (true); \ } \ else \ { \ r = SparseBoolMatrix (nr, nc, m.nnz ()); \ r.cidx (0) = static_cast<octave_idx_type> (0); \ octave_idx_type nel = 0; \ for (octave_idx_type j = 0; j < nc; j++) \ { \ for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ if ((m.data (i) != LHS_ZERO) OP (s != RHS_ZERO)) \ { \ r.ridx (nel) = m.ridx (i); \ r.data (nel++) = true; \ } \ r.cidx (j + 1) = nel; \ } \ r.maybe_compress (false); \ } \ } \ return r; \ } #define SPARSE_SMS_BOOL_OPS2(M, S, LHS_ZERO, RHS_ZERO) \ SPARSE_SMS_BOOL_OP (mx_el_and, &&, M, S, LHS_ZERO, RHS_ZERO) \ SPARSE_SMS_BOOL_OP (mx_el_or, ||, M, S, LHS_ZERO, RHS_ZERO) #define SPARSE_SMS_BOOL_OPS(M, S, ZERO) \ SPARSE_SMS_BOOL_OPS2(M, S, ZERO, ZERO) // scalar by sparse matrix operations. #define SPARSE_SSM_BIN_OP_1(R, F, OP, S, M) \ R \ F (const S& s, const M& m) \ { \ octave_idx_type nr = m.rows (); \ octave_idx_type nc = m.cols (); \ \ R r (nr, nc, (s OP 0.0)); \ \ for (octave_idx_type j = 0; j < nc; j++) \ for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ r.xelem (m.ridx (i), j) = s OP m.data (i); \ \ return r; \ } #define SPARSE_SSM_BIN_OP_2(R, F, OP, S, M) \ R \ F (const S& s, const M& m) \ { \ octave_idx_type nr = m.rows (); \ octave_idx_type nc = m.cols (); \ octave_idx_type nz = m.nnz (); \ \ R r (nr, nc, nz); \ \ for (octave_idx_type i = 0; i < nz; i++) \ { \ r.xdata (i) = s OP m.data (i); \ r.xridx (i) = m.ridx (i); \ } \ for (octave_idx_type i = 0; i < nc + 1; i++) \ r.xcidx (i) = m.cidx (i); \ \ r.maybe_compress(true); \ return r; \ } #define SPARSE_SSM_BIN_OPS(R1, R2, S, M) \ SPARSE_SSM_BIN_OP_1 (R1, operator +, +, S, M) \ SPARSE_SSM_BIN_OP_1 (R1, operator -, -, S, M) \ SPARSE_SSM_BIN_OP_2 (R2, operator *, *, S, M) \ SPARSE_SSM_BIN_OP_2 (R2, operator /, /, S, M) #define SPARSE_SSM_CMP_OP(F, OP, S, SZ, SC, M, MZ, MC) \ SparseBoolMatrix \ F (const S& s, const M& m) \ { \ octave_idx_type nr = m.rows (); \ octave_idx_type nc = m.cols (); \ SparseBoolMatrix r; \ \ if (SC (s) OP SC (MZ)) \ { \ r = SparseBoolMatrix (nr, nc, true); \ for (octave_idx_type j = 0; j < nc; j++) \ for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ if (! (SC (s) OP MC (m.data (i)))) \ r.data (m.ridx (i) + j * nr) = false; \ r.maybe_compress (true); \ } \ else \ { \ r = SparseBoolMatrix (nr, nc, m.nnz ()); \ r.cidx (0) = static_cast<octave_idx_type> (0); \ octave_idx_type nel = 0; \ for (octave_idx_type j = 0; j < nc; j++) \ { \ for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ if (SC (s) OP MC (m.data (i))) \ { \ r.ridx (nel) = m.ridx (i); \ r.data (nel++) = true; \ } \ r.cidx (j + 1) = nel; \ } \ r.maybe_compress (false); \ } \ return r; \ } #define SPARSE_SSM_CMP_OPS(S, SZ, SC, M, MZ, MC) \ SPARSE_SSM_CMP_OP (mx_el_lt, <, S, SZ, , M, MZ, ) \ SPARSE_SSM_CMP_OP (mx_el_le, <=, S, SZ, , M, MZ, ) \ SPARSE_SSM_CMP_OP (mx_el_ge, >=, S, SZ, , M, MZ, ) \ SPARSE_SSM_CMP_OP (mx_el_gt, >, S, SZ, , M, MZ, ) \ SPARSE_SSM_CMP_OP (mx_el_eq, ==, S, SZ, , M, MZ, ) \ SPARSE_SSM_CMP_OP (mx_el_ne, !=, S, SZ, , M, MZ, ) #define SPARSE_SSM_EQNE_OPS(S, SZ, SC, M, MZ, MC) \ SPARSE_SSM_CMP_OP (mx_el_eq, ==, S, SZ, , M, MZ, ) \ SPARSE_SSM_CMP_OP (mx_el_ne, !=, S, SZ, , M, MZ, ) #define SPARSE_SSM_BOOL_OP(F, OP, S, M, LHS_ZERO, RHS_ZERO) \ SparseBoolMatrix \ F (const S& s, const M& m) \ { \ octave_idx_type nr = m.rows (); \ octave_idx_type nc = m.cols (); \ SparseBoolMatrix r; \ \ if (nr > 0 && nc > 0) \ { \ if ((s != LHS_ZERO) OP RHS_ZERO) \ { \ r = SparseBoolMatrix (nr, nc, true); \ for (octave_idx_type j = 0; j < nc; j++) \ for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ if (! ((s != LHS_ZERO) OP (m.data (i) != RHS_ZERO))) \ r.data (m.ridx (i) + j * nr) = false; \ r.maybe_compress (true); \ } \ else \ { \ r = SparseBoolMatrix (nr, nc, m.nnz ()); \ r.cidx (0) = static_cast<octave_idx_type> (0); \ octave_idx_type nel = 0; \ for (octave_idx_type j = 0; j < nc; j++) \ { \ for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \ if ((s != LHS_ZERO) OP (m.data (i) != RHS_ZERO)) \ { \ r.ridx (nel) = m.ridx (i); \ r.data (nel++) = true; \ } \ r.cidx (j + 1) = nel; \ } \ r.maybe_compress (false); \ } \ } \ return r; \ } #define SPARSE_SSM_BOOL_OPS2(S, M, LHS_ZERO, RHS_ZERO) \ SPARSE_SSM_BOOL_OP (mx_el_and, &&, S, M, LHS_ZERO, RHS_ZERO) \ SPARSE_SSM_BOOL_OP (mx_el_or, ||, S, M, LHS_ZERO, RHS_ZERO) #define SPARSE_SSM_BOOL_OPS(S, M, ZERO) \ SPARSE_SSM_BOOL_OPS2(S, M, ZERO, ZERO) // sparse matrix by sparse matrix operations. #define SPARSE_SMSM_BIN_OP_1(R, F, OP, M1, M2) \ R \ F (const M1& m1, const M2& m2) \ { \ R r; \ \ octave_idx_type m1_nr = m1.rows (); \ octave_idx_type m1_nc = m1.cols (); \ \ octave_idx_type m2_nr = m2.rows (); \ octave_idx_type m2_nc = m2.cols (); \ \ if (m1_nr == 1 && m1_nc == 1) \ { \ if (m1.elem (0,0) == 0.) \ r = OP R (m2); \ else \ { \ r = R (m2_nr, m2_nc, m1.data (0) OP 0.); \ \ for (octave_idx_type j = 0 ; j < m2_nc ; j++) \ { \ octave_quit (); \ octave_idx_type idxj = j * m2_nr; \ for (octave_idx_type i = m2.cidx (j) ; i < m2.cidx (j+1) ; i++) \ { \ octave_quit (); \ r.data (idxj + m2.ridx (i)) = m1.data (0) OP m2.data (i); \ } \ } \ r.maybe_compress (); \ } \ } \ else if (m2_nr == 1 && m2_nc == 1) \ { \ if (m2.elem (0,0) == 0.) \ r = R (m1); \ else \ { \ r = R (m1_nr, m1_nc, 0. OP m2.data (0)); \ \ for (octave_idx_type j = 0 ; j < m1_nc ; j++) \ { \ octave_quit (); \ octave_idx_type idxj = j * m1_nr; \ for (octave_idx_type i = m1.cidx (j) ; i < m1.cidx (j+1) ; i++) \ { \ octave_quit (); \ r.data (idxj + m1.ridx (i)) = m1.data (i) OP m2.data (0); \ } \ } \ r.maybe_compress (); \ } \ } \ else if (m1_nr != m2_nr || m1_nc != m2_nc) \ err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ else \ { \ r = R (m1_nr, m1_nc, (m1.nnz () + m2.nnz ())); \ \ octave_idx_type jx = 0; \ r.cidx (0) = 0; \ for (octave_idx_type i = 0 ; i < m1_nc ; i++) \ { \ octave_idx_type ja = m1.cidx (i); \ octave_idx_type ja_max = m1.cidx (i+1); \ bool ja_lt_max= ja < ja_max; \ \ octave_idx_type jb = m2.cidx (i); \ octave_idx_type jb_max = m2.cidx (i+1); \ bool jb_lt_max = jb < jb_max; \ \ while (ja_lt_max || jb_lt_max) \ { \ octave_quit (); \ if ((! jb_lt_max) || \ (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \ { \ r.ridx (jx) = m1.ridx (ja); \ r.data (jx) = m1.data (ja) OP 0.; \ jx++; \ ja++; \ ja_lt_max= ja < ja_max; \ } \ else if ((! ja_lt_max) || \ (jb_lt_max && (m2.ridx (jb) < m1.ridx (ja)))) \ { \ r.ridx (jx) = m2.ridx (jb); \ r.data (jx) = 0. OP m2.data (jb); \ jx++; \ jb++; \ jb_lt_max= jb < jb_max; \ } \ else \ { \ if ((m1.data (ja) OP m2.data (jb)) != 0.) \ { \ r.data (jx) = m1.data (ja) OP m2.data (jb); \ r.ridx (jx) = m1.ridx (ja); \ jx++; \ } \ ja++; \ ja_lt_max= ja < ja_max; \ jb++; \ jb_lt_max= jb < jb_max; \ } \ } \ r.cidx (i+1) = jx; \ } \ \ r.maybe_compress (); \ } \ \ return r; \ } #define SPARSE_SMSM_BIN_OP_2(R, F, OP, M1, M2) \ R \ F (const M1& m1, const M2& m2) \ { \ R r; \ \ octave_idx_type m1_nr = m1.rows (); \ octave_idx_type m1_nc = m1.cols (); \ \ octave_idx_type m2_nr = m2.rows (); \ octave_idx_type m2_nc = m2.cols (); \ \ if (m1_nr == 1 && m1_nc == 1) \ { \ if (m1.elem (0,0) == 0.) \ r = R (m2_nr, m2_nc); \ else \ { \ r = R (m2); \ octave_idx_type m2_nnz = m2.nnz (); \ \ for (octave_idx_type i = 0 ; i < m2_nnz ; i++) \ { \ octave_quit (); \ r.data (i) = m1.data (0) OP r.data (i); \ } \ r.maybe_compress (); \ } \ } \ else if (m2_nr == 1 && m2_nc == 1) \ { \ if (m2.elem (0,0) == 0.) \ r = R (m1_nr, m1_nc); \ else \ { \ r = R (m1); \ octave_idx_type m1_nnz = m1.nnz (); \ \ for (octave_idx_type i = 0 ; i < m1_nnz ; i++) \ { \ octave_quit (); \ r.data (i) = r.data (i) OP m2.data (0); \ } \ r.maybe_compress (); \ } \ } \ else if (m1_nr != m2_nr || m1_nc != m2_nc) \ err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ else \ { \ r = R (m1_nr, m1_nc, (m1.nnz () > m2.nnz () ? m1.nnz () : m2.nnz ())); \ \ octave_idx_type jx = 0; \ r.cidx (0) = 0; \ for (octave_idx_type i = 0 ; i < m1_nc ; i++) \ { \ octave_idx_type ja = m1.cidx (i); \ octave_idx_type ja_max = m1.cidx (i+1); \ bool ja_lt_max= ja < ja_max; \ \ octave_idx_type jb = m2.cidx (i); \ octave_idx_type jb_max = m2.cidx (i+1); \ bool jb_lt_max = jb < jb_max; \ \ while (ja_lt_max || jb_lt_max) \ { \ octave_quit (); \ if ((! jb_lt_max) || \ (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \ { \ ja++; ja_lt_max= ja < ja_max; \ } \ else if ((! ja_lt_max) || \ (jb_lt_max && (m2.ridx (jb) < m1.ridx (ja)))) \ { \ jb++; jb_lt_max= jb < jb_max; \ } \ else \ { \ if ((m1.data (ja) OP m2.data (jb)) != 0.) \ { \ r.data (jx) = m1.data (ja) OP m2.data (jb); \ r.ridx (jx) = m1.ridx (ja); \ jx++; \ } \ ja++; ja_lt_max= ja < ja_max; \ jb++; jb_lt_max= jb < jb_max; \ } \ } \ r.cidx (i+1) = jx; \ } \ \ r.maybe_compress (); \ } \ \ return r; \ } #define SPARSE_SMSM_BIN_OP_3(R, F, OP, M1, M2) \ R \ F (const M1& m1, const M2& m2) \ { \ R r; \ \ octave_idx_type m1_nr = m1.rows (); \ octave_idx_type m1_nc = m1.cols (); \ \ octave_idx_type m2_nr = m2.rows (); \ octave_idx_type m2_nc = m2.cols (); \ \ if (m1_nr == 1 && m1_nc == 1) \ { \ if ((m1.elem (0,0) OP Complex ()) == Complex ()) \ { \ octave_idx_type m2_nnz = m2.nnz (); \ r = R (m2); \ for (octave_idx_type i = 0 ; i < m2_nnz ; i++) \ r.data (i) = m1.elem (0,0) OP r.data (i); \ r.maybe_compress (); \ } \ else \ { \ r = R (m2_nr, m2_nc, m1.elem (0,0) OP Complex ()); \ for (octave_idx_type j = 0 ; j < m2_nc ; j++) \ { \ octave_quit (); \ octave_idx_type idxj = j * m2_nr; \ for (octave_idx_type i = m2.cidx (j) ; i < m2.cidx (j+1) ; i++) \ { \ octave_quit (); \ r.data (idxj + m2.ridx (i)) = m1.elem (0,0) OP m2.data (i); \ } \ } \ r.maybe_compress (); \ } \ } \ else if (m2_nr == 1 && m2_nc == 1) \ { \ if ((Complex () OP m1.elem (0,0)) == Complex ()) \ { \ octave_idx_type m1_nnz = m1.nnz (); \ r = R (m1); \ for (octave_idx_type i = 0 ; i < m1_nnz ; i++) \ r.data (i) = r.data (i) OP m2.elem (0,0); \ r.maybe_compress (); \ } \ else \ { \ r = R (m1_nr, m1_nc, Complex () OP m2.elem (0,0)); \ for (octave_idx_type j = 0 ; j < m1_nc ; j++) \ { \ octave_quit (); \ octave_idx_type idxj = j * m1_nr; \ for (octave_idx_type i = m1.cidx (j) ; i < m1.cidx (j+1) ; i++) \ { \ octave_quit (); \ r.data (idxj + m1.ridx (i)) = m1.data (i) OP m2.elem (0,0); \ } \ } \ r.maybe_compress (); \ } \ } \ else if (m1_nr != m2_nr || m1_nc != m2_nc) \ err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ else \ { \ \ /* FIXME: Kludge... Always double/Complex, so Complex () */ \ r = R (m1_nr, m1_nc, (Complex () OP Complex ())); \ \ for (octave_idx_type i = 0 ; i < m1_nc ; i++) \ { \ octave_idx_type ja = m1.cidx (i); \ octave_idx_type ja_max = m1.cidx (i+1); \ bool ja_lt_max= ja < ja_max; \ \ octave_idx_type jb = m2.cidx (i); \ octave_idx_type jb_max = m2.cidx (i+1); \ bool jb_lt_max = jb < jb_max; \ \ while (ja_lt_max || jb_lt_max) \ { \ octave_quit (); \ if ((! jb_lt_max) || \ (ja_lt_max && (m1.ridx (ja) < m2.ridx (jb)))) \ { \ /* keep those kludges coming */ \ r.elem (m1.ridx (ja),i) = m1.data (ja) OP Complex (); \ ja++; \ ja_lt_max= ja < ja_max; \ } \ else if ((! ja_lt_max) || \ (jb_lt_max && (m2.ridx (jb) < m1.ridx (ja)))) \ { \ /* keep those kludges coming */ \ r.elem (m2.ridx (jb),i) = Complex () OP m2.data (jb); \ jb++; \ jb_lt_max= jb < jb_max; \ } \ else \ { \ r.elem (m1.ridx (ja),i) = m1.data (ja) OP m2.data (jb); \ ja++; \ ja_lt_max= ja < ja_max; \ jb++; \ jb_lt_max= jb < jb_max; \ } \ } \ } \ r.maybe_compress (true); \ } \ \ return r; \ } // Note that SM ./ SM needs to take into account the NaN and Inf values // implied by the division by zero. // FIXME: Are the NaNs double(NaN) or Complex(NaN,Nan) in the complex case? #define SPARSE_SMSM_BIN_OPS(R1, R2, M1, M2) \ SPARSE_SMSM_BIN_OP_1 (R1, operator +, +, M1, M2) \ SPARSE_SMSM_BIN_OP_1 (R1, operator -, -, M1, M2) \ SPARSE_SMSM_BIN_OP_2 (R2, product, *, M1, M2) \ SPARSE_SMSM_BIN_OP_3 (R2, quotient, /, M1, M2) // FIXME: this macro duplicates the bodies of the template functions // defined in the SPARSE_SSM_CMP_OP and SPARSE_SMS_CMP_OP macros. #define SPARSE_SMSM_CMP_OP(F, OP, M1, Z1, C1, M2, Z2, C2) \ SparseBoolMatrix \ F (const M1& m1, const M2& m2) \ { \ SparseBoolMatrix r; \ \ octave_idx_type m1_nr = m1.rows (); \ octave_idx_type m1_nc = m1.cols (); \ \ octave_idx_type m2_nr = m2.rows (); \ octave_idx_type m2_nc = m2.cols (); \ \ if (m1_nr == 1 && m1_nc == 1) \ { \ if (C1 (m1.elem (0,0)) OP C2 (Z2)) \ { \ r = SparseBoolMatrix (m2_nr, m2_nc, true); \ for (octave_idx_type j = 0; j < m2_nc; j++) \ for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \ if (! (C1 (m1.elem (0,0)) OP C2 (m2.data (i)))) \ r.data (m2.ridx (i) + j * m2_nr) = false; \ r.maybe_compress (true); \ } \ else \ { \ r = SparseBoolMatrix (m2_nr, m2_nc, m2.nnz ()); \ r.cidx (0) = static_cast<octave_idx_type> (0); \ octave_idx_type nel = 0; \ for (octave_idx_type j = 0; j < m2_nc; j++) \ { \ for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \ if (C1 (m1.elem (0,0)) OP C2 (m2.data (i))) \ { \ r.ridx (nel) = m2.ridx (i); \ r.data (nel++) = true; \ } \ r.cidx (j + 1) = nel; \ } \ r.maybe_compress (false); \ } \ } \ else if (m2_nr == 1 && m2_nc == 1) \ { \ if (C1 (Z1) OP C2 (m2.elem (0,0))) \ { \ r = SparseBoolMatrix (m1_nr, m1_nc, true); \ for (octave_idx_type j = 0; j < m1_nc; j++) \ for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \ if (! (C1 (m1.data (i)) OP C2 (m2.elem (0,0)))) \ r.data (m1.ridx (i) + j * m1_nr) = false; \ r.maybe_compress (true); \ } \ else \ { \ r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz ()); \ r.cidx (0) = static_cast<octave_idx_type> (0); \ octave_idx_type nel = 0; \ for (octave_idx_type j = 0; j < m1_nc; j++) \ { \ for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \ if (C1 (m1.data (i)) OP C2 (m2.elem (0,0))) \ { \ r.ridx (nel) = m1.ridx (i); \ r.data (nel++) = true; \ } \ r.cidx (j + 1) = nel; \ } \ r.maybe_compress (false); \ } \ } \ else if (m1_nr == m2_nr && m1_nc == m2_nc) \ { \ if (m1_nr != 0 || m1_nc != 0) \ { \ if (C1 (Z1) OP C2 (Z2)) \ { \ r = SparseBoolMatrix (m1_nr, m1_nc, true); \ for (octave_idx_type j = 0; j < m1_nc; j++) \ { \ octave_idx_type i1 = m1.cidx (j); \ octave_idx_type e1 = m1.cidx (j+1); \ octave_idx_type i2 = m2.cidx (j); \ octave_idx_type e2 = m2.cidx (j+1); \ while (i1 < e1 || i2 < e2) \ { \ if (i1 == e1 || (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \ { \ if (! (C1 (Z1) OP C2 (m2.data (i2)))) \ r.data (m2.ridx (i2) + j * m1_nr) = false; \ i2++; \ } \ else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \ { \ if (! (C1 (m1.data (i1)) OP C2 (Z2))) \ r.data (m1.ridx (i1) + j * m1_nr) = false; \ i1++; \ } \ else \ { \ if (! (C1 (m1.data (i1)) OP C2 (m2.data (i2)))) \ r.data (m1.ridx (i1) + j * m1_nr) = false; \ i1++; \ i2++; \ } \ } \ } \ r.maybe_compress (true); \ } \ else \ { \ r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz () + m2.nnz ()); \ r.cidx (0) = static_cast<octave_idx_type> (0); \ octave_idx_type nel = 0; \ for (octave_idx_type j = 0; j < m1_nc; j++) \ { \ octave_idx_type i1 = m1.cidx (j); \ octave_idx_type e1 = m1.cidx (j+1); \ octave_idx_type i2 = m2.cidx (j); \ octave_idx_type e2 = m2.cidx (j+1); \ while (i1 < e1 || i2 < e2) \ { \ if (i1 == e1 || (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \ { \ if (C1 (Z1) OP C2 (m2.data (i2))) \ { \ r.ridx (nel) = m2.ridx (i2); \ r.data (nel++) = true; \ } \ i2++; \ } \ else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \ { \ if (C1 (m1.data (i1)) OP C2 (Z2)) \ { \ r.ridx (nel) = m1.ridx (i1); \ r.data (nel++) = true; \ } \ i1++; \ } \ else \ { \ if (C1 (m1.data (i1)) OP C2 (m2.data (i2))) \ { \ r.ridx (nel) = m1.ridx (i1); \ r.data (nel++) = true; \ } \ i1++; \ i2++; \ } \ } \ r.cidx (j + 1) = nel; \ } \ r.maybe_compress (false); \ } \ } \ } \ else \ { \ if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ } \ return r; \ } #define SPARSE_SMSM_CMP_OPS(M1, Z1, C1, M2, Z2, C2) \ SPARSE_SMSM_CMP_OP (mx_el_lt, <, M1, Z1, , M2, Z2, ) \ SPARSE_SMSM_CMP_OP (mx_el_le, <=, M1, Z1, , M2, Z2, ) \ SPARSE_SMSM_CMP_OP (mx_el_ge, >=, M1, Z1, , M2, Z2, ) \ SPARSE_SMSM_CMP_OP (mx_el_gt, >, M1, Z1, , M2, Z2, ) \ SPARSE_SMSM_CMP_OP (mx_el_eq, ==, M1, Z1, , M2, Z2, ) \ SPARSE_SMSM_CMP_OP (mx_el_ne, !=, M1, Z1, , M2, Z2, ) #define SPARSE_SMSM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2) \ SPARSE_SMSM_CMP_OP (mx_el_eq, ==, M1, Z1, , M2, Z2, ) \ SPARSE_SMSM_CMP_OP (mx_el_ne, !=, M1, Z1, , M2, Z2, ) // FIXME: this macro duplicates the bodies of the template functions // defined in the SPARSE_SSM_BOOL_OP and SPARSE_SMS_BOOL_OP macros. #define SPARSE_SMSM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \ SparseBoolMatrix \ F (const M1& m1, const M2& m2) \ { \ SparseBoolMatrix r; \ \ octave_idx_type m1_nr = m1.rows (); \ octave_idx_type m1_nc = m1.cols (); \ \ octave_idx_type m2_nr = m2.rows (); \ octave_idx_type m2_nc = m2.cols (); \ \ if (m1_nr == 1 && m1_nc == 1) \ { \ if (m2_nr > 0 && m2_nc > 0) \ { \ if ((m1.elem (0,0) != LHS_ZERO) OP RHS_ZERO) \ { \ r = SparseBoolMatrix (m2_nr, m2_nc, true); \ for (octave_idx_type j = 0; j < m2_nc; j++) \ for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \ if (! ((m1.elem (0,0) != LHS_ZERO) OP (m2.data (i) != RHS_ZERO))) \ r.data (m2.ridx (i) + j * m2_nr) = false; \ r.maybe_compress (true); \ } \ else \ { \ r = SparseBoolMatrix (m2_nr, m2_nc, m2.nnz ()); \ r.cidx (0) = static_cast<octave_idx_type> (0); \ octave_idx_type nel = 0; \ for (octave_idx_type j = 0; j < m2_nc; j++) \ { \ for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \ if ((m1.elem (0,0) != LHS_ZERO) OP (m2.data (i) != RHS_ZERO)) \ { \ r.ridx (nel) = m2.ridx (i); \ r.data (nel++) = true; \ } \ r.cidx (j + 1) = nel; \ } \ r.maybe_compress (false); \ } \ } \ } \ else if (m2_nr == 1 && m2_nc == 1) \ { \ if (m1_nr > 0 && m1_nc > 0) \ { \ if (LHS_ZERO OP (m2.elem (0,0) != RHS_ZERO)) \ { \ r = SparseBoolMatrix (m1_nr, m1_nc, true); \ for (octave_idx_type j = 0; j < m1_nc; j++) \ for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \ if (! ((m1.data (i) != LHS_ZERO) OP (m2.elem (0,0) != RHS_ZERO))) \ r.data (m1.ridx (i) + j * m1_nr) = false; \ r.maybe_compress (true); \ } \ else \ { \ r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz ()); \ r.cidx (0) = static_cast<octave_idx_type> (0); \ octave_idx_type nel = 0; \ for (octave_idx_type j = 0; j < m1_nc; j++) \ { \ for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \ if ((m1.data (i) != LHS_ZERO) OP (m2.elem (0,0) != RHS_ZERO)) \ { \ r.ridx (nel) = m1.ridx (i); \ r.data (nel++) = true; \ } \ r.cidx (j + 1) = nel; \ } \ r.maybe_compress (false); \ } \ } \ } \ else if (m1_nr == m2_nr && m1_nc == m2_nc) \ { \ if (m1_nr != 0 || m1_nc != 0) \ { \ r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz () + m2.nnz ()); \ r.cidx (0) = static_cast<octave_idx_type> (0); \ octave_idx_type nel = 0; \ for (octave_idx_type j = 0; j < m1_nc; j++) \ { \ octave_idx_type i1 = m1.cidx (j); \ octave_idx_type e1 = m1.cidx (j+1); \ octave_idx_type i2 = m2.cidx (j); \ octave_idx_type e2 = m2.cidx (j+1); \ while (i1 < e1 || i2 < e2) \ { \ if (i1 == e1 || (i2 < e2 && m1.ridx (i1) > m2.ridx (i2))) \ { \ if (LHS_ZERO OP m2.data (i2) != RHS_ZERO) \ { \ r.ridx (nel) = m2.ridx (i2); \ r.data (nel++) = true; \ } \ i2++; \ } \ else if (i2 == e2 || m1.ridx (i1) < m2.ridx (i2)) \ { \ if (m1.data (i1) != LHS_ZERO OP RHS_ZERO) \ { \ r.ridx (nel) = m1.ridx (i1); \ r.data (nel++) = true; \ } \ i1++; \ } \ else \ { \ if (m1.data (i1) != LHS_ZERO OP m2.data (i2) != RHS_ZERO) \ { \ r.ridx (nel) = m1.ridx (i1); \ r.data (nel++) = true; \ } \ i1++; \ i2++; \ } \ } \ r.cidx (j + 1) = nel; \ } \ r.maybe_compress (false); \ } \ } \ else \ { \ if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ } \ return r; \ } #define SPARSE_SMSM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \ SPARSE_SMSM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \ SPARSE_SMSM_BOOL_OP (mx_el_or, ||, M1, M2, LHS_ZERO, RHS_ZERO) #define SPARSE_SMSM_BOOL_OPS(M1, M2, ZERO) \ SPARSE_SMSM_BOOL_OPS2(M1, M2, ZERO, ZERO) // matrix by sparse matrix operations. #define SPARSE_MSM_BIN_OP_1(R, F, OP, M1, M2) \ R \ F (const M1& m1, const M2& m2) \ { \ R r; \ \ octave_idx_type m1_nr = m1.rows (); \ octave_idx_type m1_nc = m1.cols (); \ \ octave_idx_type m2_nr = m2.rows (); \ octave_idx_type m2_nc = m2.cols (); \ \ if (m2_nr == 1 && m2_nc == 1) \ r = R (m1 OP m2.elem (0,0)); \ else if (m1_nr != m2_nr || m1_nc != m2_nc) \ err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ else \ { \ r = R (F (m1, m2.matrix_value ())); \ } \ return r; \ } #define SPARSE_MSM_BIN_OP_2(R, F, OP, M1, M2) \ R \ F (const M1& m1, const M2& m2) \ { \ R r; \ \ octave_idx_type m1_nr = m1.rows (); \ octave_idx_type m1_nc = m1.cols (); \ \ octave_idx_type m2_nr = m2.rows (); \ octave_idx_type m2_nc = m2.cols (); \ \ if (m2_nr == 1 && m2_nc == 1) \ r = R (m1 OP m2.elem (0,0)); \ else if (m1_nr != m2_nr || m1_nc != m2_nc) \ err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ else \ { \ if (do_mx_check (m1, mx_inline_all_finite<M1::element_type>)) \ { \ /* Sparsity pattern is preserved. */ \ octave_idx_type m2_nz = m2.nnz (); \ r = R (m2_nr, m2_nc, m2_nz); \ for (octave_idx_type j = 0, k = 0; j < m2_nc; j++) \ { \ octave_quit (); \ for (octave_idx_type i = m2.cidx (j); i < m2.cidx (j+1); i++) \ { \ octave_idx_type mri = m2.ridx (i); \ R::element_type x = m1(mri, j) OP m2.data (i); \ if (x != 0.0) \ { \ r.xdata (k) = x; \ r.xridx (k) = m2.ridx (i); \ k++; \ } \ } \ r.xcidx (j+1) = k; \ } \ r.maybe_compress (false); \ return r; \ } \ else \ r = R (F (m1, m2.matrix_value ())); \ } \ \ return r; \ } // FIXME: Pass a specific ZERO value #define SPARSE_MSM_BIN_OPS(R1, R2, M1, M2) \ SPARSE_MSM_BIN_OP_1 (R1, operator +, +, M1, M2) \ SPARSE_MSM_BIN_OP_1 (R1, operator -, -, M1, M2) \ SPARSE_MSM_BIN_OP_2 (R2, product, *, M1, M2) \ SPARSE_MSM_BIN_OP_1 (R2, quotient, /, M1, M2) #define SPARSE_MSM_CMP_OP(F, OP, M1, C1, M2, C2) \ SparseBoolMatrix \ F (const M1& m1, const M2& m2) \ { \ SparseBoolMatrix r; \ \ octave_idx_type m1_nr = m1.rows (); \ octave_idx_type m1_nc = m1.cols (); \ \ octave_idx_type m2_nr = m2.rows (); \ octave_idx_type m2_nc = m2.cols (); \ \ if (m2_nr == 1 && m2_nc == 1) \ r = SparseBoolMatrix (F (m1, m2.elem (0,0))); \ else if (m1_nr == m2_nr && m1_nc == m2_nc) \ { \ if (m1_nr != 0 || m1_nc != 0) \ { \ /* Count num of nonzero elements */ \ octave_idx_type nel = 0; \ for (octave_idx_type j = 0; j < m1_nc; j++) \ for (octave_idx_type i = 0; i < m1_nr; i++) \ if (C1 (m1.elem (i, j)) OP C2 (m2.elem (i, j))) \ nel++; \ \ r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ \ octave_idx_type ii = 0; \ r.cidx (0) = 0; \ for (octave_idx_type j = 0; j < m1_nc; j++) \ { \ for (octave_idx_type i = 0; i < m1_nr; i++) \ { \ bool el = C1 (m1.elem (i, j)) OP C2 (m2.elem (i, j)); \ if (el) \ { \ r.data (ii) = el; \ r.ridx (ii++) = i; \ } \ } \ r.cidx (j+1) = ii; \ } \ } \ } \ else \ { \ if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ } \ return r; \ } #define SPARSE_MSM_CMP_OPS(M1, Z1, C1, M2, Z2, C2) \ SPARSE_MSM_CMP_OP (mx_el_lt, <, M1, , M2, ) \ SPARSE_MSM_CMP_OP (mx_el_le, <=, M1, , M2, ) \ SPARSE_MSM_CMP_OP (mx_el_ge, >=, M1, , M2, ) \ SPARSE_MSM_CMP_OP (mx_el_gt, >, M1, , M2, ) \ SPARSE_MSM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \ SPARSE_MSM_CMP_OP (mx_el_ne, !=, M1, , M2, ) #define SPARSE_MSM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2) \ SPARSE_MSM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \ SPARSE_MSM_CMP_OP (mx_el_ne, !=, M1, , M2, ) #define SPARSE_MSM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \ SparseBoolMatrix \ F (const M1& m1, const M2& m2) \ { \ SparseBoolMatrix r; \ \ octave_idx_type m1_nr = m1.rows (); \ octave_idx_type m1_nc = m1.cols (); \ \ octave_idx_type m2_nr = m2.rows (); \ octave_idx_type m2_nc = m2.cols (); \ \ if (m2_nr == 1 && m2_nc == 1) \ r = SparseBoolMatrix (F (m1, m2.elem (0,0))); \ else if (m1_nr == m2_nr && m1_nc == m2_nc) \ { \ if (m1_nr != 0 || m1_nc != 0) \ { \ /* Count num of nonzero elements */ \ octave_idx_type nel = 0; \ for (octave_idx_type j = 0; j < m1_nc; j++) \ for (octave_idx_type i = 0; i < m1_nr; i++) \ if ((m1.elem (i, j) != LHS_ZERO) \ OP (m2.elem (i, j) != RHS_ZERO)) \ nel++; \ \ r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ \ octave_idx_type ii = 0; \ r.cidx (0) = 0; \ for (octave_idx_type j = 0; j < m1_nc; j++) \ { \ for (octave_idx_type i = 0; i < m1_nr; i++) \ { \ bool el = (m1.elem (i, j) != LHS_ZERO) \ OP (m2.elem (i, j) != RHS_ZERO); \ if (el) \ { \ r.data (ii) = el; \ r.ridx (ii++) = i; \ } \ } \ r.cidx (j+1) = ii; \ } \ } \ } \ else \ { \ if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ } \ return r; \ } #define SPARSE_MSM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \ SPARSE_MSM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \ SPARSE_MSM_BOOL_OP (mx_el_or, ||, M1, M2, LHS_ZERO, RHS_ZERO) #define SPARSE_MSM_BOOL_OPS(M1, M2, ZERO) \ SPARSE_MSM_BOOL_OPS2(M1, M2, ZERO, ZERO) // sparse matrix by matrix operations. #define SPARSE_SMM_BIN_OP_1(R, F, OP, M1, M2) \ R \ F (const M1& m1, const M2& m2) \ { \ R r; \ \ octave_idx_type m1_nr = m1.rows (); \ octave_idx_type m1_nc = m1.cols (); \ \ octave_idx_type m2_nr = m2.rows (); \ octave_idx_type m2_nc = m2.cols (); \ \ if (m1_nr == 1 && m1_nc == 1) \ r = R (m1.elem (0,0) OP m2); \ else if (m1_nr != m2_nr || m1_nc != m2_nc) \ err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ else \ { \ r = R (m1.matrix_value () OP m2); \ } \ return r; \ } // sm .* m preserves sparsity if m contains no Infs nor Nans. #define SPARSE_SMM_BIN_OP_2_CHECK_product(ET) \ do_mx_check (m2, mx_inline_all_finite<ET>) // sm ./ m preserves sparsity if m contains no NaNs or zeros. #define SPARSE_SMM_BIN_OP_2_CHECK_quotient(ET) \ ! do_mx_check (m2, mx_inline_any_nan<ET>) && m2.nnz () == m2.numel () #define SPARSE_SMM_BIN_OP_2(R, F, OP, M1, M2) \ R \ F (const M1& m1, const M2& m2) \ { \ R r; \ \ octave_idx_type m1_nr = m1.rows (); \ octave_idx_type m1_nc = m1.cols (); \ \ octave_idx_type m2_nr = m2.rows (); \ octave_idx_type m2_nc = m2.cols (); \ \ if (m1_nr == 1 && m1_nc == 1) \ r = R (m1.elem (0,0) OP m2); \ else if (m1_nr != m2_nr || m1_nc != m2_nc) \ err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ else \ { \ if (SPARSE_SMM_BIN_OP_2_CHECK_ ## F(M2::element_type)) \ { \ /* Sparsity pattern is preserved. */ \ octave_idx_type m1_nz = m1.nnz (); \ r = R (m1_nr, m1_nc, m1_nz); \ for (octave_idx_type j = 0, k = 0; j < m1_nc; j++) \ { \ octave_quit (); \ for (octave_idx_type i = m1.cidx (j); i < m1.cidx (j+1); i++) \ { \ octave_idx_type mri = m1.ridx (i); \ R::element_type x = m1.data (i) OP m2 (mri, j); \ if (x != 0.0) \ { \ r.xdata (k) = x; \ r.xridx (k) = m1.ridx (i); \ k++; \ } \ } \ r.xcidx (j+1) = k; \ } \ r.maybe_compress (false); \ return r; \ } \ else \ r = R (F (m1.matrix_value (), m2)); \ } \ \ return r; \ } #define SPARSE_SMM_BIN_OPS(R1, R2, M1, M2) \ SPARSE_SMM_BIN_OP_1 (R1, operator +, +, M1, M2) \ SPARSE_SMM_BIN_OP_1 (R1, operator -, -, M1, M2) \ SPARSE_SMM_BIN_OP_2 (R2, product, *, M1, M2) \ SPARSE_SMM_BIN_OP_2 (R2, quotient, /, M1, M2) #define SPARSE_SMM_CMP_OP(F, OP, M1, C1, M2, C2) \ SparseBoolMatrix \ F (const M1& m1, const M2& m2) \ { \ SparseBoolMatrix r; \ \ octave_idx_type m1_nr = m1.rows (); \ octave_idx_type m1_nc = m1.cols (); \ \ octave_idx_type m2_nr = m2.rows (); \ octave_idx_type m2_nc = m2.cols (); \ \ if (m1_nr == 1 && m1_nc == 1) \ r = SparseBoolMatrix (F (m1.elem (0,0), m2)); \ else if (m1_nr == m2_nr && m1_nc == m2_nc) \ { \ if (m1_nr != 0 || m1_nc != 0) \ { \ /* Count num of nonzero elements */ \ octave_idx_type nel = 0; \ for (octave_idx_type j = 0; j < m1_nc; j++) \ for (octave_idx_type i = 0; i < m1_nr; i++) \ if (C1 (m1.elem (i, j)) OP C2 (m2.elem (i, j))) \ nel++; \ \ r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ \ octave_idx_type ii = 0; \ r.cidx (0) = 0; \ for (octave_idx_type j = 0; j < m1_nc; j++) \ { \ for (octave_idx_type i = 0; i < m1_nr; i++) \ { \ bool el = C1 (m1.elem (i, j)) OP C2 (m2.elem (i, j)); \ if (el) \ { \ r.data (ii) = el; \ r.ridx (ii++) = i; \ } \ } \ r.cidx (j+1) = ii; \ } \ } \ } \ else \ { \ if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ } \ return r; \ } #define SPARSE_SMM_CMP_OPS(M1, Z1, C1, M2, Z2, C2) \ SPARSE_SMM_CMP_OP (mx_el_lt, <, M1, , M2, ) \ SPARSE_SMM_CMP_OP (mx_el_le, <=, M1, , M2, ) \ SPARSE_SMM_CMP_OP (mx_el_ge, >=, M1, , M2, ) \ SPARSE_SMM_CMP_OP (mx_el_gt, >, M1, , M2, ) \ SPARSE_SMM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \ SPARSE_SMM_CMP_OP (mx_el_ne, !=, M1, , M2, ) #define SPARSE_SMM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2) \ SPARSE_SMM_CMP_OP (mx_el_eq, ==, M1, , M2, ) \ SPARSE_SMM_CMP_OP (mx_el_ne, !=, M1, , M2, ) #define SPARSE_SMM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \ SparseBoolMatrix \ F (const M1& m1, const M2& m2) \ { \ SparseBoolMatrix r; \ \ octave_idx_type m1_nr = m1.rows (); \ octave_idx_type m1_nc = m1.cols (); \ \ octave_idx_type m2_nr = m2.rows (); \ octave_idx_type m2_nc = m2.cols (); \ \ if (m1_nr == 1 && m1_nc == 1) \ r = SparseBoolMatrix (F (m1.elem (0,0), m2)); \ else if (m1_nr == m2_nr && m1_nc == m2_nc) \ { \ if (m1_nr != 0 || m1_nc != 0) \ { \ /* Count num of nonzero elements */ \ octave_idx_type nel = 0; \ for (octave_idx_type j = 0; j < m1_nc; j++) \ for (octave_idx_type i = 0; i < m1_nr; i++) \ if ((m1.elem (i, j) != LHS_ZERO) \ OP (m2.elem (i, j) != RHS_ZERO)) \ nel++; \ \ r = SparseBoolMatrix (m1_nr, m1_nc, nel); \ \ octave_idx_type ii = 0; \ r.cidx (0) = 0; \ for (octave_idx_type j = 0; j < m1_nc; j++) \ { \ for (octave_idx_type i = 0; i < m1_nr; i++) \ { \ bool el = (m1.elem (i, j) != LHS_ZERO) \ OP (m2.elem (i, j) != RHS_ZERO); \ if (el) \ { \ r.data (ii) = el; \ r.ridx (ii++) = i; \ } \ } \ r.cidx (j+1) = ii; \ } \ } \ } \ else \ { \ if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \ err_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \ } \ return r; \ } #define SPARSE_SMM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \ SPARSE_SMM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \ SPARSE_SMM_BOOL_OP (mx_el_or, ||, M1, M2, LHS_ZERO, RHS_ZERO) #define SPARSE_SMM_BOOL_OPS(M1, M2, ZERO) \ SPARSE_SMM_BOOL_OPS2(M1, M2, ZERO, ZERO) // Avoid some code duplication. Maybe we should use templates. #define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN) \ \ octave_idx_type nr = rows (); \ octave_idx_type nc = cols (); \ \ RET_TYPE retval; \ \ if (nr > 0 && nc > 0) \ { \ if ((nr == 1 && dim == -1) || dim == 1) \ /* Ugly!! Is there a better way? */ \ retval = transpose (). FCN (0) .transpose (); \ else \ { \ octave_idx_type nel = 0; \ for (octave_idx_type i = 0; i < nc; i++) \ { \ ELT_TYPE t = ELT_TYPE (); \ for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \ { \ t += data (j); \ if (t != ELT_TYPE ()) \ { \ if (j == cidx (i+1) - 1) \ nel += nr - ridx (j); \ else \ nel += ridx (j+1) - ridx (j); \ } \ } \ } \ retval = RET_TYPE (nr, nc, nel); \ retval.cidx (0) = 0; \ octave_idx_type ii = 0; \ for (octave_idx_type i = 0; i < nc; i++) \ { \ ELT_TYPE t = ELT_TYPE (); \ for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \ { \ t += data (j); \ if (t != ELT_TYPE ()) \ { \ if (j == cidx (i+1) - 1) \ { \ for (octave_idx_type k = ridx (j); k < nr; k++) \ { \ retval.data (ii) = t; \ retval.ridx (ii++) = k; \ } \ } \ else \ { \ for (octave_idx_type k = ridx (j); k < ridx (j+1); k++) \ { \ retval.data (ii) = t; \ retval.ridx (ii++) = k; \ } \ } \ } \ } \ retval.cidx (i+1) = ii; \ } \ } \ } \ else \ retval = RET_TYPE (nr,nc); \ \ return retval #define SPARSE_CUMPROD(RET_TYPE, ELT_TYPE, FCN) \ \ octave_idx_type nr = rows (); \ octave_idx_type nc = cols (); \ \ RET_TYPE retval; \ \ if (nr > 0 && nc > 0) \ { \ if ((nr == 1 && dim == -1) || dim == 1) \ /* Ugly!! Is there a better way? */ \ retval = transpose (). FCN (0) .transpose (); \ else \ { \ octave_idx_type nel = 0; \ for (octave_idx_type i = 0; i < nc; i++) \ { \ octave_idx_type jj = 0; \ for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \ { \ if (jj == ridx (j)) \ { \ nel++; \ jj++; \ } \ else \ break; \ } \ } \ retval = RET_TYPE (nr, nc, nel); \ retval.cidx (0) = 0; \ octave_idx_type ii = 0; \ for (octave_idx_type i = 0; i < nc; i++) \ { \ ELT_TYPE t = ELT_TYPE (1.); \ octave_idx_type jj = 0; \ for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \ { \ if (jj == ridx (j)) \ { \ t *= data (j); \ retval.data (ii) = t; \ retval.ridx (ii++) = jj++; \ } \ else \ break; \ } \ retval.cidx (i+1) = ii; \ } \ } \ } \ else \ retval = RET_TYPE (nr,nc); \ \ return retval #define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, \ INIT_VAL, MT_RESULT) \ \ octave_idx_type nr = rows (); \ octave_idx_type nc = cols (); \ \ RET_TYPE retval; \ \ if (nr > 0 && nc > 0) \ { \ if ((nr == 1 && dim == -1) || dim == 1) \ { \ /* Define j here to allow fancy definition for prod method */ \ octave_idx_type j = 0; \ OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nr); \ \ for (octave_idx_type i = 0; i < nr; i++) \ tmp[i] = INIT_VAL; \ for (j = 0; j < nc; j++) \ { \ for (octave_idx_type i = cidx (j); i < cidx (j + 1); i++) \ { \ ROW_EXPR; \ } \ } \ octave_idx_type nel = 0; \ for (octave_idx_type i = 0; i < nr; i++) \ if (tmp[i] != EL_TYPE ()) \ nel++ ; \ retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), nel); \ retval.cidx (0) = 0; \ retval.cidx (1) = nel; \ nel = 0; \ for (octave_idx_type i = 0; i < nr; i++) \ if (tmp[i] != EL_TYPE ()) \ { \ retval.data (nel) = tmp[i]; \ retval.ridx (nel++) = i; \ } \ } \ else \ { \ OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nc); \ \ for (octave_idx_type j = 0; j < nc; j++) \ { \ tmp[j] = INIT_VAL; \ for (octave_idx_type i = cidx (j); i < cidx (j + 1); i++) \ { \ COL_EXPR; \ } \ } \ octave_idx_type nel = 0; \ for (octave_idx_type i = 0; i < nc; i++) \ if (tmp[i] != EL_TYPE ()) \ nel++ ; \ retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, nel); \ retval.cidx (0) = 0; \ nel = 0; \ for (octave_idx_type i = 0; i < nc; i++) \ if (tmp[i] != EL_TYPE ()) \ { \ retval.data (nel) = tmp[i]; \ retval.ridx (nel++) = 0; \ retval.cidx (i+1) = retval.cidx (i) + 1; \ } \ else \ retval.cidx (i+1) = retval.cidx (i); \ } \ } \ else if (nc == 0 && (nr == 0 || (nr == 1 && dim == -1))) \ { \ if (MT_RESULT) \ { \ retval = RET_TYPE (static_cast<octave_idx_type> (1), \ static_cast<octave_idx_type> (1), \ static_cast<octave_idx_type> (1)); \ retval.cidx (0) = 0; \ retval.cidx (1) = 1; \ retval.ridx (0) = 0; \ retval.data (0) = MT_RESULT; \ } \ else \ retval = RET_TYPE (static_cast<octave_idx_type> (1), \ static_cast<octave_idx_type> (1), \ static_cast<octave_idx_type> (0)); \ } \ else if (nr == 0 && (dim == 0 || dim == -1)) \ { \ if (MT_RESULT) \ { \ retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, nc); \ retval.cidx (0) = 0; \ for (octave_idx_type i = 0; i < nc ; i++) \ { \ retval.ridx (i) = 0; \ retval.cidx (i+1) = i+1; \ retval.data (i) = MT_RESULT; \ } \ } \ else \ retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, \ static_cast<octave_idx_type> (0)); \ } \ else if (nc == 0 && dim == 1) \ { \ if (MT_RESULT) \ { \ retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), nr); \ retval.cidx (0) = 0; \ retval.cidx (1) = nr; \ for (octave_idx_type i = 0; i < nr; i++) \ { \ retval.ridx (i) = i; \ retval.data (i) = MT_RESULT; \ } \ } \ else \ retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), \ static_cast<octave_idx_type> (0)); \ } \ else \ retval.resize (nr > 0, nc > 0); \ \ return retval #define SPARSE_REDUCTION_OP_ROW_EXPR(OP) \ tmp[ridx (i)] OP data (i) #define SPARSE_REDUCTION_OP_COL_EXPR(OP) \ tmp[j] OP data (i) #define SPARSE_REDUCTION_OP(RET_TYPE, EL_TYPE, OP, INIT_VAL, MT_RESULT) \ SPARSE_BASE_REDUCTION_OP (RET_TYPE, EL_TYPE, \ SPARSE_REDUCTION_OP_ROW_EXPR (OP), \ SPARSE_REDUCTION_OP_COL_EXPR (OP), \ INIT_VAL, MT_RESULT) // Don't break from this loop if the test succeeds because // we are looping over the rows and not the columns in the inner loop. #define SPARSE_ANY_ALL_OP_ROW_CODE(TEST_OP, TEST_TRUE_VAL) \ if (data (i) TEST_OP 0.0) \ tmp[ridx (i)] = TEST_TRUE_VAL; #define SPARSE_ANY_ALL_OP_COL_CODE(TEST_OP, TEST_TRUE_VAL) \ if (data (i) TEST_OP 0.0) \ { \ tmp[j] = TEST_TRUE_VAL; \ break; \ } #define SPARSE_ANY_ALL_OP(DIM, INIT_VAL, MT_RESULT, TEST_OP, TEST_TRUE_VAL) \ SPARSE_BASE_REDUCTION_OP (SparseBoolMatrix, char, \ SPARSE_ANY_ALL_OP_ROW_CODE (TEST_OP, TEST_TRUE_VAL), \ SPARSE_ANY_ALL_OP_COL_CODE (TEST_OP, TEST_TRUE_VAL), \ INIT_VAL, MT_RESULT) #define SPARSE_ALL_OP(DIM) \ if ((rows () == 1 && dim == -1) || dim == 1) \ return transpose (). all (0). transpose (); \ else \ { \ SPARSE_ANY_ALL_OP (DIM, (cidx (j+1) - cidx (j) < nr ? false : true), \ true, ==, false); \ } #define SPARSE_ANY_OP(DIM) SPARSE_ANY_ALL_OP (DIM, false, false, !=, true) #define SPARSE_SPARSE_MUL(RET_TYPE, RET_EL_TYPE, EL_TYPE) \ octave_idx_type nr = m.rows (); \ octave_idx_type nc = m.cols (); \ \ octave_idx_type a_nr = a.rows (); \ octave_idx_type a_nc = a.cols (); \ \ if (nr == 1 && nc == 1) \ { \ RET_EL_TYPE s = m.elem (0,0); \ octave_idx_type nz = a.nnz (); \ RET_TYPE r (a_nr, a_nc, nz); \ \ for (octave_idx_type i = 0; i < nz; i++) \ { \ octave_quit (); \ r.data (i) = s * a.data (i); \ r.ridx (i) = a.ridx (i); \ } \ for (octave_idx_type i = 0; i < a_nc + 1; i++) \ { \ octave_quit (); \ r.cidx (i) = a.cidx (i); \ } \ \ r.maybe_compress (true); \ return r; \ } \ else if (a_nr == 1 && a_nc == 1) \ { \ RET_EL_TYPE s = a.elem (0,0); \ octave_idx_type nz = m.nnz (); \ RET_TYPE r (nr, nc, nz); \ \ for (octave_idx_type i = 0; i < nz; i++) \ { \ octave_quit (); \ r.data (i) = m.data (i) * s; \ r.ridx (i) = m.ridx (i); \ } \ for (octave_idx_type i = 0; i < nc + 1; i++) \ { \ octave_quit (); \ r.cidx (i) = m.cidx (i); \ } \ \ r.maybe_compress (true); \ return r; \ } \ else if (nc != a_nr) \ err_nonconformant ("operator *", nr, nc, a_nr, a_nc); \ else \ { \ OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr); \ RET_TYPE retval (nr, a_nc, static_cast<octave_idx_type> (0)); \ for (octave_idx_type i = 0; i < nr; i++) \ w[i] = 0; \ retval.xcidx (0) = 0; \ \ octave_idx_type nel = 0; \ \ for (octave_idx_type i = 0; i < a_nc; i++) \ { \ for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \ { \ octave_idx_type col = a.ridx (j); \ for (octave_idx_type k = m.cidx (col) ; k < m.cidx (col+1); k++) \ { \ if (w[m.ridx (k)] < i + 1) \ { \ w[m.ridx (k)] = i + 1; \ nel++; \ } \ octave_quit (); \ } \ } \ retval.xcidx (i+1) = nel; \ } \ \ if (nel == 0) \ return RET_TYPE (nr, a_nc); \ else \ { \ for (octave_idx_type i = 0; i < nr; i++) \ w[i] = 0; \ \ OCTAVE_LOCAL_BUFFER (RET_EL_TYPE, Xcol, nr); \ \ retval.change_capacity (nel); \ /* The optimal break-point as estimated from simulations */ \ /* Note that Mergesort is O(nz log(nz)) while searching all */ \ /* values is O(nr), where nz here is nonzero per row of */ \ /* length nr. The test itself was then derived from the */ \ /* simulation with random square matrices and the observation */ \ /* of the number of nonzero elements in the output matrix */ \ /* it was found that the breakpoints were */ \ /* nr: 500 1000 2000 5000 10000 */ \ /* nz: 6 25 97 585 2202 */ \ /* The below is a simplication of the 'polyfit'-ed parameters */ \ /* to these breakpoints */ \ octave_idx_type n_per_col = (a_nc > 43000 ? 43000 : \ (a_nc * a_nc) / 43000); \ octave_idx_type ii = 0; \ octave_idx_type *ri = retval.xridx (); \ octave_sort<octave_idx_type> sort; \ \ for (octave_idx_type i = 0; i < a_nc ; i++) \ { \ if (retval.xcidx (i+1) - retval.xcidx (i) > n_per_col) \ { \ for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \ { \ octave_idx_type col = a.ridx (j); \ EL_TYPE tmpval = a.data (j); \ for (octave_idx_type k = m.cidx (col) ; \ k < m.cidx (col+1); k++) \ { \ octave_quit (); \ octave_idx_type row = m.ridx (k); \ if (w[row] < i + 1) \ { \ w[row] = i + 1; \ Xcol[row] = tmpval * m.data (k); \ } \ else \ Xcol[row] += tmpval * m.data (k); \ } \ } \ for (octave_idx_type k = 0; k < nr; k++) \ if (w[k] == i + 1) \ { \ retval.xdata (ii) = Xcol[k]; \ retval.xridx (ii++) = k; \ } \ } \ else \ { \ for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \ { \ octave_idx_type col = a.ridx (j); \ EL_TYPE tmpval = a.data (j); \ for (octave_idx_type k = m.cidx (col) ; \ k < m.cidx (col+1); k++) \ { \ octave_quit (); \ octave_idx_type row = m.ridx (k); \ if (w[row] < i + 1) \ { \ w[row] = i + 1; \ retval.xridx (ii++) = row;\ Xcol[row] = tmpval * m.data (k); \ } \ else \ Xcol[row] += tmpval * m.data (k); \ } \ } \ sort.sort (ri + retval.xcidx (i), ii - retval.xcidx (i)); \ for (octave_idx_type k = retval.xcidx (i); k < ii; k++) \ retval.xdata (k) = Xcol[retval.xridx (k)]; \ } \ } \ retval.maybe_compress (true);\ return retval; \ } \ } #define SPARSE_FULL_MUL(RET_TYPE, EL_TYPE, ZERO) \ octave_idx_type nr = m.rows (); \ octave_idx_type nc = m.cols (); \ \ octave_idx_type a_nr = a.rows (); \ octave_idx_type a_nc = a.cols (); \ \ if (nr == 1 && nc == 1) \ { \ RET_TYPE retval = m.elem (0,0) * a; \ return retval; \ } \ else if (nc != a_nr) \ err_nonconformant ("operator *", nr, nc, a_nr, a_nc); \ else \ { \ RET_TYPE retval (nr, a_nc, ZERO); \ \ for (octave_idx_type i = 0; i < a_nc ; i++) \ { \ for (octave_idx_type j = 0; j < a_nr; j++) \ { \ octave_quit (); \ \ EL_TYPE tmpval = a.elem (j,i); \ for (octave_idx_type k = m.cidx (j) ; k < m.cidx (j+1); k++) \ retval.elem (m.ridx (k),i) += tmpval * m.data (k); \ } \ } \ return retval; \ } #define SPARSE_FULL_TRANS_MUL(RET_TYPE, EL_TYPE, ZERO, CONJ_OP) \ octave_idx_type nr = m.rows (); \ octave_idx_type nc = m.cols (); \ \ octave_idx_type a_nr = a.rows (); \ octave_idx_type a_nc = a.cols (); \ \ if (nr == 1 && nc == 1) \ { \ RET_TYPE retval = CONJ_OP (m.elem (0,0)) * a; \ return retval; \ } \ else if (nr != a_nr) \ err_nonconformant ("operator *", nc, nr, a_nr, a_nc); \ else \ { \ RET_TYPE retval (nc, a_nc); \ \ for (octave_idx_type i = 0; i < a_nc ; i++) \ { \ for (octave_idx_type j = 0; j < nc; j++) \ { \ octave_quit (); \ \ EL_TYPE acc = ZERO; \ for (octave_idx_type k = m.cidx (j) ; k < m.cidx (j+1); k++) \ acc += a.elem (m.ridx (k),i) * CONJ_OP (m.data (k)); \ retval.xelem (j,i) = acc; \ } \ } \ return retval; \ } #define FULL_SPARSE_MUL(RET_TYPE, EL_TYPE, ZERO) \ octave_idx_type nr = m.rows (); \ octave_idx_type nc = m.cols (); \ \ octave_idx_type a_nr = a.rows (); \ octave_idx_type a_nc = a.cols (); \ \ if (a_nr == 1 && a_nc == 1) \ { \ RET_TYPE retval = m * a.elem (0,0); \ return retval; \ } \ else if (nc != a_nr) \ err_nonconformant ("operator *", nr, nc, a_nr, a_nc); \ else \ { \ RET_TYPE retval (nr, a_nc, ZERO); \ \ for (octave_idx_type i = 0; i < a_nc ; i++) \ { \ octave_quit (); \ for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \ { \ octave_idx_type col = a.ridx (j); \ EL_TYPE tmpval = a.data (j); \ \ for (octave_idx_type k = 0 ; k < nr; k++) \ retval.xelem (k,i) += tmpval * m.elem (k,col); \ } \ } \ return retval; \ } #define FULL_SPARSE_MUL_TRANS(RET_TYPE, EL_TYPE, ZERO, CONJ_OP) \ octave_idx_type nr = m.rows (); \ octave_idx_type nc = m.cols (); \ \ octave_idx_type a_nr = a.rows (); \ octave_idx_type a_nc = a.cols (); \ \ if (a_nr == 1 && a_nc == 1) \ { \ RET_TYPE retval = m * CONJ_OP (a.elem (0,0)); \ return retval; \ } \ else if (nc != a_nc) \ err_nonconformant ("operator *", nr, nc, a_nc, a_nr); \ else \ { \ RET_TYPE retval (nr, a_nr, ZERO); \ \ for (octave_idx_type i = 0; i < a_nc ; i++) \ { \ octave_quit (); \ for (octave_idx_type j = a.cidx (i); j < a.cidx (i+1); j++) \ { \ octave_idx_type col = a.ridx (j); \ EL_TYPE tmpval = CONJ_OP (a.data (j)); \ for (octave_idx_type k = 0 ; k < nr; k++) \ retval.xelem (k,col) += tmpval * m.elem (k,i); \ } \ } \ return retval; \ } #endif