Sparse-op-defs.h

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 2004-2012 David Bateman
00004 Copyright (C) 1998-2004 Andy Adler
00005 Copyright (C) 2008 Jaroslav Hajek
00006 
00007 This file is part of Octave.
00008 
00009 Octave is free software; you can redistribute it and/or modify it
00010 under the terms of the GNU General Public License as published by the
00011 Free Software Foundation; either version 3 of the License, or (at your
00012 option) any later version.
00013 
00014 Octave is distributed in the hope that it will be useful, but WITHOUT
00015 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00016 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00017 for more details.
00018 
00019 You should have received a copy of the GNU General Public License
00020 along with Octave; see the file COPYING.  If not, see
00021 <http://www.gnu.org/licenses/>.
00022 
00023 */
00024 
00025 #if !defined (octave_sparse_op_defs_h)
00026 #define octave_sparse_op_defs_h 1
00027 
00028 #include "Array-util.h"
00029 #include "mx-ops.h"
00030 #include "oct-locbuf.h"
00031 #include "mx-inlines.cc"
00032 
00033 #define SPARSE_BIN_OP_DECL(R, OP, X, Y, API) \
00034   extern API R OP (const X&, const Y&)
00035 
00036 #define SPARSE_CMP_OP_DECL(OP, X, Y, API) \
00037   extern API SparseBoolMatrix OP (const X&, const Y&)
00038 
00039 #define SPARSE_BOOL_OP_DECL(OP, X, Y, API) \
00040   extern API SparseBoolMatrix OP (const X&, const Y&)
00041 
00042 // matrix by scalar operations.
00043 
00044 #define SPARSE_SMS_BIN_OP_DECLS(R1, R2, M, S, API)  \
00045   SPARSE_BIN_OP_DECL (R1, operator +, M, S, API); \
00046   SPARSE_BIN_OP_DECL (R1, operator -, M, S, API); \
00047   SPARSE_BIN_OP_DECL (R2, operator *, M, S, API); \
00048   SPARSE_BIN_OP_DECL (R2, operator /, M, S, API);
00049 
00050 #define SPARSE_SMS_BIN_OP_1(R, F, OP, M, S)     \
00051   R \
00052   F (const M& m, const S& s) \
00053   { \
00054     octave_idx_type nr = m.rows (); \
00055     octave_idx_type nc = m.cols (); \
00056  \
00057     R r (nr, nc, (0.0 OP s)); \
00058  \
00059     for (octave_idx_type j = 0; j < nc; j++) \
00060       for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
00061         r.xelem (m.ridx (i), j) = m.data (i) OP s; \
00062     return r; \
00063   }
00064 
00065 #define SPARSE_SMS_BIN_OP_2(R, F, OP, M, S)     \
00066   R \
00067   F (const M& m, const S& s) \
00068   { \
00069     octave_idx_type nr = m.rows (); \
00070     octave_idx_type nc = m.cols (); \
00071     octave_idx_type nz = m.nnz (); \
00072  \
00073     R r (nr, nc, nz); \
00074  \
00075     for (octave_idx_type i = 0; i < nz; i++) \
00076       { \
00077         r.xdata(i) = m.data(i) OP s; \
00078         r.xridx(i) = m.ridx(i); \
00079       } \
00080     for (octave_idx_type i = 0; i < nc + 1; i++) \
00081       r.xcidx(i) = m.cidx(i); \
00082     \
00083     r.maybe_compress (true); \
00084     return r; \
00085   }
00086 
00087 #define SPARSE_SMS_BIN_OPS(R1, R2, M, S) \
00088   SPARSE_SMS_BIN_OP_1 (R1, operator +, +, M, S) \
00089   SPARSE_SMS_BIN_OP_1 (R1, operator -, -, M, S) \
00090   SPARSE_SMS_BIN_OP_2 (R2, operator *, *, M, S) \
00091   SPARSE_SMS_BIN_OP_2 (R2, operator /, /, M, S)
00092 
00093 #define SPARSE_SMS_CMP_OP_DECLS(M, S, API) \
00094   SPARSE_CMP_OP_DECL (mx_el_lt, M, S, API); \
00095   SPARSE_CMP_OP_DECL (mx_el_le, M, S, API); \
00096   SPARSE_CMP_OP_DECL (mx_el_ge, M, S, API); \
00097   SPARSE_CMP_OP_DECL (mx_el_gt, M, S, API); \
00098   SPARSE_CMP_OP_DECL (mx_el_eq, M, S, API); \
00099   SPARSE_CMP_OP_DECL (mx_el_ne, M, S, API);
00100 
00101 #define SPARSE_SMS_EQNE_OP_DECLS(M, S, API) \
00102   SPARSE_CMP_OP_DECL (mx_el_eq, M, S, API); \
00103   SPARSE_CMP_OP_DECL (mx_el_ne, M, S, API);
00104 
00105 #define SPARSE_SMS_CMP_OP(F, OP, M, MZ, MC, S, SZ, SC)  \
00106   SparseBoolMatrix \
00107   F (const M& m, const S& s) \
00108   { \
00109     octave_idx_type nr = m.rows (); \
00110     octave_idx_type nc = m.cols (); \
00111     SparseBoolMatrix r; \
00112     \
00113     if (MC (MZ) OP SC (s)) \
00114       { \
00115         r = SparseBoolMatrix (nr, nc, true); \
00116         for (octave_idx_type j = 0; j < nc; j++) \
00117           for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) \
00118             if (! (MC (m.data (i)) OP SC (s))) \
00119               r.data (m.ridx (i) + j * nr) = false; \
00120         r.maybe_compress (true); \
00121       } \
00122     else \
00123       { \
00124         r = SparseBoolMatrix (nr, nc, m.nnz ()); \
00125         r.cidx (0) = static_cast<octave_idx_type> (0); \
00126         octave_idx_type nel = 0; \
00127         for (octave_idx_type j = 0; j < nc; j++) \
00128           { \
00129             for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) \
00130               if (MC (m.data (i)) OP SC (s)) \
00131                 { \
00132                   r.ridx (nel) = m.ridx (i); \
00133                   r.data (nel++) = true; \
00134                 } \
00135             r.cidx (j + 1) = nel; \
00136           } \
00137         r.maybe_compress (false); \
00138       } \
00139     return r; \
00140   }
00141 
00142 #define SPARSE_SMS_CMP_OPS(M, MZ, CM, S, SZ, CS)        \
00143   SPARSE_SMS_CMP_OP (mx_el_lt, <,  M, MZ,   , S, SZ,   )        \
00144   SPARSE_SMS_CMP_OP (mx_el_le, <=, M, MZ,   , S, SZ,   )        \
00145   SPARSE_SMS_CMP_OP (mx_el_ge, >=, M, MZ,   , S, SZ,   )        \
00146   SPARSE_SMS_CMP_OP (mx_el_gt, >,  M, MZ,   , S, SZ,   )        \
00147   SPARSE_SMS_CMP_OP (mx_el_eq, ==, M, MZ,   , S, SZ,   )        \
00148   SPARSE_SMS_CMP_OP (mx_el_ne, !=, M, MZ,   , S, SZ,   )
00149 
00150 #define SPARSE_SMS_EQNE_OPS(M, MZ, CM, S, SZ, CS)       \
00151   SPARSE_SMS_CMP_OP (mx_el_eq, ==, M, MZ,   , S, SZ,   )        \
00152   SPARSE_SMS_CMP_OP (mx_el_ne, !=, M, MZ,   , S, SZ,   )
00153 
00154 #define SPARSE_SMS_BOOL_OP_DECLS(M, S, API) \
00155   SPARSE_BOOL_OP_DECL (mx_el_and, M, S, API); \
00156   SPARSE_BOOL_OP_DECL (mx_el_or,  M, S, API);
00157 
00158 #define SPARSE_SMS_BOOL_OP(F, OP, M, S, LHS_ZERO, RHS_ZERO) \
00159   SparseBoolMatrix \
00160   F (const M& m, const S& s) \
00161   { \
00162     octave_idx_type nr = m.rows (); \
00163     octave_idx_type nc = m.cols (); \
00164     SparseBoolMatrix r; \
00165     \
00166     if (nr > 0 && nc > 0) \
00167       { \
00168         if (LHS_ZERO OP (s != RHS_ZERO)) \
00169           { \
00170             r = SparseBoolMatrix (nr, nc, true); \
00171             for (octave_idx_type j = 0; j < nc; j++) \
00172               for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) \
00173                 if (! ((m.data(i) != LHS_ZERO) OP (s != RHS_ZERO))) \
00174                   r.data (m.ridx (i) + j * nr) = false; \
00175             r.maybe_compress (true); \
00176           } \
00177         else \
00178           { \
00179             r = SparseBoolMatrix (nr, nc, m.nnz ()); \
00180             r.cidx (0) = static_cast<octave_idx_type> (0); \
00181             octave_idx_type nel = 0; \
00182             for (octave_idx_type j = 0; j < nc; j++) \
00183               { \
00184                 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) \
00185                   if ((m.data(i) != LHS_ZERO) OP (s != RHS_ZERO)) \
00186                     { \
00187                       r.ridx (nel) = m.ridx (i); \
00188                       r.data (nel++) = true; \
00189                     } \
00190                 r.cidx (j + 1) = nel; \
00191               } \
00192             r.maybe_compress (false); \
00193           } \
00194       } \
00195     return r; \
00196   }
00197 
00198 #define SPARSE_SMS_BOOL_OPS2(M, S, LHS_ZERO, RHS_ZERO) \
00199   SPARSE_SMS_BOOL_OP (mx_el_and, &&, M, S, LHS_ZERO, RHS_ZERO) \
00200   SPARSE_SMS_BOOL_OP (mx_el_or,  ||, M, S, LHS_ZERO, RHS_ZERO)
00201 
00202 #define SPARSE_SMS_BOOL_OPS(M, S, ZERO) \
00203   SPARSE_SMS_BOOL_OPS2(M, S, ZERO, ZERO)
00204 
00205 #define SPARSE_SMS_OP_DECLS(R1, R2, M, S, API) \
00206   SPARSE_SMS_BIN_OP_DECLS (R1, R2, M, S, API)    \
00207   SPARSE_SMS_CMP_OP_DECLS (M, S, API) \
00208   SPARSE_SMS_BOOL_OP_DECLS (M, S, API)
00209 
00210 // scalar by matrix operations.
00211 
00212 #define SPARSE_SSM_BIN_OP_DECLS(R1, R2, S, M, API)    \
00213   SPARSE_BIN_OP_DECL (R1, operator +, S, M, API); \
00214   SPARSE_BIN_OP_DECL (R1, operator -, S, M, API); \
00215   SPARSE_BIN_OP_DECL (R2, operator *, S, M, API); \
00216   SPARSE_BIN_OP_DECL (R2, operator /, S, M, API);
00217 
00218 #define SPARSE_SSM_BIN_OP_1(R, F, OP, S, M) \
00219   R \
00220   F (const S& s, const M& m) \
00221   { \
00222     octave_idx_type nr = m.rows (); \
00223     octave_idx_type nc = m.cols (); \
00224  \
00225     R r (nr, nc, (s OP 0.0)); \
00226  \
00227     for (octave_idx_type j = 0; j < nc; j++) \
00228       for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++) \
00229         r.xelem (m.ridx (i), j) = s OP m.data (i); \
00230  \
00231     return r; \
00232   }
00233 
00234 #define SPARSE_SSM_BIN_OP_2(R, F, OP, S, M) \
00235   R \
00236   F (const S& s, const M& m) \
00237   { \
00238     octave_idx_type nr = m.rows (); \
00239     octave_idx_type nc = m.cols (); \
00240     octave_idx_type nz = m.nnz (); \
00241  \
00242     R r (nr, nc, nz); \
00243  \
00244     for (octave_idx_type i = 0; i < nz; i++) \
00245       { \
00246         r.xdata(i) = s OP m.data(i); \
00247         r.xridx(i) = m.ridx(i); \
00248       } \
00249     for (octave_idx_type i = 0; i < nc + 1; i++) \
00250       r.xcidx(i) = m.cidx(i); \
00251  \
00252     r.maybe_compress(true); \
00253     return r; \
00254   }
00255 
00256 #define SPARSE_SSM_BIN_OPS(R1, R2, S, M) \
00257   SPARSE_SSM_BIN_OP_1 (R1, operator +, +, S, M) \
00258   SPARSE_SSM_BIN_OP_1 (R1, operator -, -, S, M) \
00259   SPARSE_SSM_BIN_OP_2 (R2, operator *, *, S, M) \
00260   SPARSE_SSM_BIN_OP_2 (R2, operator /, /, S, M)
00261 
00262 #define SPARSE_SSM_CMP_OP_DECLS(S, M, API) \
00263   SPARSE_CMP_OP_DECL (mx_el_lt, S, M, API); \
00264   SPARSE_CMP_OP_DECL (mx_el_le, S, M, API); \
00265   SPARSE_CMP_OP_DECL (mx_el_ge, S, M, API); \
00266   SPARSE_CMP_OP_DECL (mx_el_gt, S, M, API); \
00267   SPARSE_CMP_OP_DECL (mx_el_eq, S, M, API); \
00268   SPARSE_CMP_OP_DECL (mx_el_ne, S, M, API);
00269 
00270 #define SPARSE_SSM_EQNE_OP_DECLS(S, M, API) \
00271   SPARSE_CMP_OP_DECL (mx_el_eq, S, M, API); \
00272   SPARSE_CMP_OP_DECL (mx_el_ne, S, M, API);
00273 
00274 #define SPARSE_SSM_CMP_OP(F, OP, S, SZ, SC, M, MZ, MC)  \
00275   SparseBoolMatrix \
00276   F (const S& s, const M& m) \
00277   { \
00278     octave_idx_type nr = m.rows (); \
00279     octave_idx_type nc = m.cols (); \
00280     SparseBoolMatrix r; \
00281     \
00282     if (SC (s) OP SC (MZ)) \
00283       { \
00284         r = SparseBoolMatrix (nr, nc, true); \
00285         for (octave_idx_type j = 0; j < nc; j++) \
00286           for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) \
00287             if (! (SC (s) OP MC (m.data (i)))) \
00288               r.data (m.ridx (i) + j * nr) = false; \
00289         r.maybe_compress (true); \
00290       } \
00291     else \
00292       { \
00293         r = SparseBoolMatrix (nr, nc, m.nnz ()); \
00294         r.cidx (0) = static_cast<octave_idx_type> (0); \
00295         octave_idx_type nel = 0; \
00296         for (octave_idx_type j = 0; j < nc; j++) \
00297           { \
00298             for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) \
00299               if (SC (s) OP MC (m.data (i))) \
00300                 { \
00301                   r.ridx (nel) = m.ridx (i); \
00302                   r.data (nel++) = true; \
00303                 } \
00304             r.cidx (j + 1) = nel; \
00305           } \
00306         r.maybe_compress (false); \
00307       } \
00308     return r; \
00309   }
00310 
00311 #define SPARSE_SSM_CMP_OPS(S, SZ, SC, M, MZ, MC)        \
00312   SPARSE_SSM_CMP_OP (mx_el_lt, <,  S, SZ,   , M, MZ,   )        \
00313   SPARSE_SSM_CMP_OP (mx_el_le, <=, S, SZ,   , M, MZ,   )        \
00314   SPARSE_SSM_CMP_OP (mx_el_ge, >=, S, SZ,   , M, MZ,   )        \
00315   SPARSE_SSM_CMP_OP (mx_el_gt, >,  S, SZ,   , M, MZ,   )        \
00316   SPARSE_SSM_CMP_OP (mx_el_eq, ==, S, SZ,   , M, MZ,   )        \
00317   SPARSE_SSM_CMP_OP (mx_el_ne, !=, S, SZ,   , M, MZ,   )
00318 
00319 #define SPARSE_SSM_EQNE_OPS(S, SZ, SC, M, MZ, MC)       \
00320   SPARSE_SSM_CMP_OP (mx_el_eq, ==, S, SZ,   , M, MZ,   )        \
00321   SPARSE_SSM_CMP_OP (mx_el_ne, !=, S, SZ,   , M, MZ,   )
00322 
00323 #define SPARSE_SSM_BOOL_OP_DECLS(S, M, API) \
00324   SPARSE_BOOL_OP_DECL (mx_el_and, S, M, API); \
00325   SPARSE_BOOL_OP_DECL (mx_el_or,  S, M, API); \
00326 
00327 #define SPARSE_SSM_BOOL_OP(F, OP, S, M, LHS_ZERO, RHS_ZERO) \
00328   SparseBoolMatrix \
00329   F (const S& s, const M& m) \
00330   { \
00331     octave_idx_type nr = m.rows (); \
00332     octave_idx_type nc = m.cols (); \
00333     SparseBoolMatrix r; \
00334     \
00335     if (nr > 0 && nc > 0) \
00336       { \
00337         if ((s != LHS_ZERO) OP RHS_ZERO) \
00338           { \
00339             r = SparseBoolMatrix (nr, nc, true); \
00340             for (octave_idx_type j = 0; j < nc; j++) \
00341               for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) \
00342                 if (! ((s != LHS_ZERO) OP (m.data(i) != RHS_ZERO))) \
00343                   r.data (m.ridx (i) + j * nr) = false; \
00344             r.maybe_compress (true); \
00345           } \
00346         else \
00347           { \
00348             r = SparseBoolMatrix (nr, nc, m.nnz ()); \
00349             r.cidx (0) = static_cast<octave_idx_type> (0); \
00350             octave_idx_type nel = 0; \
00351             for (octave_idx_type j = 0; j < nc; j++) \
00352               { \
00353                 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) \
00354                   if ((s != LHS_ZERO) OP (m.data(i) != RHS_ZERO)) \
00355                     { \
00356                       r.ridx (nel) = m.ridx (i); \
00357                       r.data (nel++) = true; \
00358                     } \
00359                 r.cidx (j + 1) = nel; \
00360               } \
00361             r.maybe_compress (false); \
00362           } \
00363       } \
00364     return r; \
00365   }
00366 
00367 #define SPARSE_SSM_BOOL_OPS2(S, M, LHS_ZERO, RHS_ZERO) \
00368   SPARSE_SSM_BOOL_OP (mx_el_and, &&, S, M, LHS_ZERO, RHS_ZERO) \
00369   SPARSE_SSM_BOOL_OP (mx_el_or,  ||, S, M, LHS_ZERO, RHS_ZERO)
00370 
00371 #define SPARSE_SSM_BOOL_OPS(S, M, ZERO) \
00372   SPARSE_SSM_BOOL_OPS2(S, M, ZERO, ZERO)
00373 
00374 #define SPARSE_SSM_OP_DECLS(R1, R2, S, M, API) \
00375   SPARSE_SSM_BIN_OP_DECLS (R1, R2, S, M, API)    \
00376   SPARSE_SSM_CMP_OP_DECLS (S, M, API) \
00377   SPARSE_SSM_BOOL_OP_DECLS (S, M, API) \
00378 
00379 // matrix by matrix operations.
00380 
00381 #define SPARSE_SMSM_BIN_OP_DECLS(R1, R2, M1, M2, API)   \
00382   SPARSE_BIN_OP_DECL (R1, operator +, M1, M2, API); \
00383   SPARSE_BIN_OP_DECL (R1, operator -, M1, M2, API); \
00384   SPARSE_BIN_OP_DECL (R2, product,    M1, M2, API); \
00385   SPARSE_BIN_OP_DECL (R2, quotient,   M1, M2, API);
00386 
00387 #define SPARSE_SMSM_BIN_OP_1(R, F, OP, M1, M2)  \
00388   R \
00389   F (const M1& m1, const M2& m2) \
00390   { \
00391     R r; \
00392  \
00393     octave_idx_type m1_nr = m1.rows (); \
00394     octave_idx_type m1_nc = m1.cols (); \
00395  \
00396     octave_idx_type m2_nr = m2.rows (); \
00397     octave_idx_type m2_nc = m2.cols (); \
00398  \
00399     if (m1_nr == 1 && m1_nc == 1) \
00400       { \
00401         if (m1.elem(0,0) == 0.) \
00402           r = OP R (m2); \
00403         else \
00404           { \
00405             r = R (m2_nr, m2_nc, m1.data(0) OP 0.); \
00406             \
00407             for (octave_idx_type j = 0 ; j < m2_nc ; j++) \
00408               { \
00409                 octave_quit (); \
00410                 octave_idx_type idxj = j * m2_nr; \
00411                 for (octave_idx_type i = m2.cidx(j) ; i < m2.cidx(j+1) ; i++) \
00412                   { \
00413                     octave_quit (); \
00414                     r.data(idxj + m2.ridx(i)) = m1.data(0) OP m2.data(i); \
00415                   } \
00416               } \
00417             r.maybe_compress (); \
00418           } \
00419       } \
00420     else if (m2_nr == 1 && m2_nc == 1) \
00421       { \
00422         if (m2.elem(0,0) == 0.) \
00423           r = R (m1); \
00424         else \
00425           { \
00426             r = R (m1_nr, m1_nc, 0. OP m2.data(0)); \
00427             \
00428             for (octave_idx_type j = 0 ; j < m1_nc ; j++) \
00429               { \
00430                 octave_quit (); \
00431                 octave_idx_type idxj = j * m1_nr; \
00432                 for (octave_idx_type i = m1.cidx(j) ; i < m1.cidx(j+1) ; i++) \
00433                   { \
00434                     octave_quit (); \
00435                     r.data(idxj + m1.ridx(i)) = m1.data(i) OP m2.data(0); \
00436                   } \
00437               } \
00438             r.maybe_compress (); \
00439           } \
00440       } \
00441     else if (m1_nr != m2_nr || m1_nc != m2_nc) \
00442       gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
00443     else \
00444       { \
00445         r = R (m1_nr, m1_nc, (m1.nnz () + m2.nnz ())); \
00446         \
00447         octave_idx_type jx = 0; \
00448         r.cidx (0) = 0; \
00449         for (octave_idx_type i = 0 ; i < m1_nc ; i++) \
00450           { \
00451             octave_idx_type  ja = m1.cidx(i); \
00452             octave_idx_type  ja_max = m1.cidx(i+1); \
00453             bool ja_lt_max= ja < ja_max; \
00454             \
00455             octave_idx_type  jb = m2.cidx(i); \
00456             octave_idx_type  jb_max = m2.cidx(i+1); \
00457             bool jb_lt_max = jb < jb_max; \
00458             \
00459             while (ja_lt_max || jb_lt_max ) \
00460               { \
00461                 octave_quit (); \
00462                 if ((! jb_lt_max) || \
00463                       (ja_lt_max && (m1.ridx(ja) < m2.ridx(jb)))) \
00464                   { \
00465                     r.ridx(jx) = m1.ridx(ja); \
00466                     r.data(jx) = m1.data(ja) OP 0.; \
00467                     jx++; \
00468                     ja++; \
00469                     ja_lt_max= ja < ja_max; \
00470                   } \
00471                 else if (( !ja_lt_max ) || \
00472                      (jb_lt_max && (m2.ridx(jb) < m1.ridx(ja)) ) ) \
00473                   { \
00474                     r.ridx(jx) = m2.ridx(jb); \
00475                     r.data(jx) = 0. OP m2.data(jb); \
00476                     jx++; \
00477                     jb++; \
00478                     jb_lt_max= jb < jb_max; \
00479                   } \
00480                 else \
00481                   { \
00482                      if ((m1.data(ja) OP m2.data(jb)) != 0.) \
00483                        { \
00484                           r.data(jx) = m1.data(ja) OP m2.data(jb); \
00485                           r.ridx(jx) = m1.ridx(ja); \
00486                           jx++; \
00487                        } \
00488                      ja++; \
00489                      ja_lt_max= ja < ja_max; \
00490                      jb++; \
00491                      jb_lt_max= jb < jb_max; \
00492                   } \
00493               } \
00494             r.cidx(i+1) = jx; \
00495           } \
00496         \
00497         r.maybe_compress (); \
00498       } \
00499  \
00500     return r; \
00501   }
00502 
00503 #define SPARSE_SMSM_BIN_OP_2(R, F, OP, M1, M2)  \
00504   R \
00505   F (const M1& m1, const M2& m2) \
00506   { \
00507     R r; \
00508  \
00509     octave_idx_type m1_nr = m1.rows (); \
00510     octave_idx_type m1_nc = m1.cols (); \
00511  \
00512     octave_idx_type m2_nr = m2.rows (); \
00513     octave_idx_type m2_nc = m2.cols (); \
00514  \
00515     if (m1_nr == 1 && m1_nc == 1) \
00516       { \
00517         if (m1.elem(0,0) == 0.) \
00518           r = R (m2_nr, m2_nc); \
00519         else \
00520           { \
00521             r = R (m2); \
00522             octave_idx_type m2_nnz = m2.nnz(); \
00523             \
00524             for (octave_idx_type i = 0 ; i < m2_nnz ; i++) \
00525               { \
00526                 octave_quit (); \
00527                 r.data (i) = m1.data(0) OP r.data(i); \
00528               } \
00529             r.maybe_compress (); \
00530           } \
00531       } \
00532     else if (m2_nr == 1 && m2_nc == 1) \
00533       { \
00534         if (m2.elem(0,0) == 0.) \
00535           r = R (m1_nr, m1_nc); \
00536         else \
00537           { \
00538             r = R (m1); \
00539             octave_idx_type m1_nnz = m1.nnz(); \
00540             \
00541             for (octave_idx_type i = 0 ; i < m1_nnz ; i++) \
00542               { \
00543                 octave_quit (); \
00544                 r.data (i) = r.data(i) OP m2.data(0); \
00545               } \
00546             r.maybe_compress (); \
00547           } \
00548       } \
00549     else if (m1_nr != m2_nr || m1_nc != m2_nc) \
00550       gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
00551     else \
00552       { \
00553         r = R (m1_nr, m1_nc, (m1.nnz () > m2.nnz () ? m1.nnz () : m2.nnz ())); \
00554         \
00555         octave_idx_type jx = 0; \
00556         r.cidx (0) = 0; \
00557         for (octave_idx_type i = 0 ; i < m1_nc ; i++) \
00558           { \
00559             octave_idx_type  ja = m1.cidx(i); \
00560             octave_idx_type  ja_max = m1.cidx(i+1); \
00561             bool ja_lt_max= ja < ja_max; \
00562             \
00563             octave_idx_type  jb = m2.cidx(i); \
00564             octave_idx_type  jb_max = m2.cidx(i+1); \
00565             bool jb_lt_max = jb < jb_max; \
00566             \
00567             while (ja_lt_max || jb_lt_max ) \
00568               { \
00569                 octave_quit (); \
00570                 if ((! jb_lt_max) || \
00571                       (ja_lt_max && (m1.ridx(ja) < m2.ridx(jb)))) \
00572                   { \
00573                      ja++; ja_lt_max= ja < ja_max; \
00574                   } \
00575                 else if (( !ja_lt_max ) || \
00576                      (jb_lt_max && (m2.ridx(jb) < m1.ridx(ja)) ) ) \
00577                   { \
00578                      jb++; jb_lt_max= jb < jb_max; \
00579                   } \
00580                 else \
00581                   { \
00582                      if ((m1.data(ja) OP m2.data(jb)) != 0.) \
00583                        { \
00584                           r.data(jx) = m1.data(ja) OP m2.data(jb); \
00585                           r.ridx(jx) = m1.ridx(ja); \
00586                           jx++; \
00587                        } \
00588                      ja++; ja_lt_max= ja < ja_max; \
00589                      jb++; jb_lt_max= jb < jb_max; \
00590                   } \
00591               } \
00592             r.cidx(i+1) = jx; \
00593           } \
00594         \
00595         r.maybe_compress (); \
00596       } \
00597  \
00598     return r; \
00599   }
00600 
00601 #define SPARSE_SMSM_BIN_OP_3(R, F, OP, M1, M2)  \
00602   R \
00603   F (const M1& m1, const M2& m2) \
00604   { \
00605     R r; \
00606  \
00607     octave_idx_type m1_nr = m1.rows (); \
00608     octave_idx_type m1_nc = m1.cols (); \
00609  \
00610     octave_idx_type m2_nr = m2.rows (); \
00611     octave_idx_type m2_nc = m2.cols (); \
00612  \
00613     if (m1_nr == 1 && m1_nc == 1) \
00614       { \
00615         if ((m1.elem (0,0) OP Complex()) == Complex()) \
00616           { \
00617             octave_idx_type m2_nnz = m2.nnz(); \
00618             r = R (m2); \
00619             for (octave_idx_type i = 0 ; i < m2_nnz ; i++) \
00620               r.data (i) = m1.elem(0,0) OP r.data(i); \
00621             r.maybe_compress (); \
00622           } \
00623         else \
00624           { \
00625             r = R (m2_nr, m2_nc, m1.elem(0,0) OP Complex ()); \
00626             for (octave_idx_type j = 0 ; j < m2_nc ; j++) \
00627               { \
00628                 octave_quit (); \
00629                 octave_idx_type idxj = j * m2_nr; \
00630                 for (octave_idx_type i = m2.cidx(j) ; i < m2.cidx(j+1) ; i++) \
00631                   { \
00632                     octave_quit (); \
00633                     r.data(idxj + m2.ridx(i)) = m1.elem(0,0) OP m2.data(i); \
00634                   } \
00635               } \
00636             r.maybe_compress (); \
00637           } \
00638       } \
00639     else if (m2_nr == 1 && m2_nc == 1) \
00640       { \
00641         if ((Complex() OP m1.elem (0,0)) == Complex()) \
00642           { \
00643             octave_idx_type m1_nnz = m1.nnz(); \
00644             r = R (m1); \
00645             for (octave_idx_type i = 0 ; i < m1_nnz ; i++) \
00646               r.data (i) = r.data(i) OP m2.elem(0,0); \
00647             r.maybe_compress (); \
00648           } \
00649         else \
00650           { \
00651             r = R (m1_nr, m1_nc, Complex() OP m2.elem(0,0)); \
00652             for (octave_idx_type j = 0 ; j < m1_nc ; j++) \
00653               { \
00654                 octave_quit (); \
00655                 octave_idx_type idxj = j * m1_nr; \
00656                 for (octave_idx_type i = m1.cidx(j) ; i < m1.cidx(j+1) ; i++) \
00657                   { \
00658                     octave_quit (); \
00659                     r.data(idxj + m1.ridx(i)) = m1.data(i) OP m2.elem(0,0); \
00660                   } \
00661               } \
00662             r.maybe_compress (); \
00663           } \
00664       } \
00665     else if (m1_nr != m2_nr || m1_nc != m2_nc) \
00666       gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
00667     else \
00668       { \
00669  \
00670         /* FIXME Kludge... Always double/Complex, so Complex () */ \
00671         r = R (m1_nr, m1_nc, (Complex () OP Complex ())); \
00672         \
00673         for (octave_idx_type i = 0 ; i < m1_nc ; i++) \
00674           { \
00675             octave_idx_type  ja = m1.cidx(i); \
00676             octave_idx_type  ja_max = m1.cidx(i+1); \
00677             bool ja_lt_max= ja < ja_max; \
00678             \
00679             octave_idx_type  jb = m2.cidx(i); \
00680             octave_idx_type  jb_max = m2.cidx(i+1); \
00681             bool jb_lt_max = jb < jb_max; \
00682             \
00683             while (ja_lt_max || jb_lt_max ) \
00684               { \
00685                 octave_quit (); \
00686                 if ((! jb_lt_max) || \
00687                       (ja_lt_max && (m1.ridx(ja) < m2.ridx(jb)))) \
00688                   { \
00689                     /* keep those kludges coming */ \
00690                     r.elem(m1.ridx(ja),i) = m1.data(ja) OP Complex (); \
00691                     ja++; \
00692                     ja_lt_max= ja < ja_max; \
00693                   } \
00694                 else if (( !ja_lt_max ) || \
00695                      (jb_lt_max && (m2.ridx(jb) < m1.ridx(ja)) ) ) \
00696                   { \
00697                     /* keep those kludges coming */ \
00698                     r.elem(m2.ridx(jb),i) = Complex () OP m2.data(jb);  \
00699                     jb++; \
00700                     jb_lt_max= jb < jb_max; \
00701                   } \
00702                 else \
00703                   { \
00704                     r.elem(m1.ridx(ja),i) = m1.data(ja) OP m2.data(jb); \
00705                     ja++; \
00706                     ja_lt_max= ja < ja_max; \
00707                     jb++; \
00708                     jb_lt_max= jb < jb_max; \
00709                   } \
00710               } \
00711           } \
00712         r.maybe_compress (true); \
00713       } \
00714  \
00715     return r; \
00716   }
00717 
00718 // Note that SM ./ SM needs to take into account the NaN and Inf values
00719 // implied by the division by zero.
00720 // FIXME Are the NaNs double(NaN) or Complex(NaN,Nan) in the complex
00721 // case?
00722 #define SPARSE_SMSM_BIN_OPS(R1, R2, M1, M2)  \
00723   SPARSE_SMSM_BIN_OP_1 (R1, operator +,  +, M1, M2) \
00724   SPARSE_SMSM_BIN_OP_1 (R1, operator -,  -, M1, M2) \
00725   SPARSE_SMSM_BIN_OP_2 (R2, product,     *, M1, M2) \
00726   SPARSE_SMSM_BIN_OP_3 (R2, quotient,    /, M1, M2)
00727 
00728 #define SPARSE_SMSM_CMP_OP_DECLS(M1, M2, API) \
00729   SPARSE_CMP_OP_DECL (mx_el_lt, M1, M2, API); \
00730   SPARSE_CMP_OP_DECL (mx_el_le, M1, M2, API); \
00731   SPARSE_CMP_OP_DECL (mx_el_ge, M1, M2, API); \
00732   SPARSE_CMP_OP_DECL (mx_el_gt, M1, M2, API); \
00733   SPARSE_CMP_OP_DECL (mx_el_eq, M1, M2, API); \
00734   SPARSE_CMP_OP_DECL (mx_el_ne, M1, M2, API);
00735 
00736 #define SPARSE_SMSM_EQNE_OP_DECLS(M1, M2, API) \
00737   SPARSE_CMP_OP_DECL (mx_el_eq, M1, M2, API); \
00738   SPARSE_CMP_OP_DECL (mx_el_ne, M1, M2, API);
00739 
00740 // FIXME -- this macro duplicatest the bodies of the template
00741 // functions defined in the SPARSE_SSM_CMP_OP and SPARSE_SMS_CMP_OP
00742 // macros.
00743 
00744 #define SPARSE_SMSM_CMP_OP(F, OP, M1, Z1, C1, M2, Z2, C2)       \
00745   SparseBoolMatrix \
00746   F (const M1& m1, const M2& m2) \
00747   { \
00748     SparseBoolMatrix r; \
00749     \
00750     octave_idx_type m1_nr = m1.rows (); \
00751     octave_idx_type m1_nc = m1.cols (); \
00752     \
00753     octave_idx_type m2_nr = m2.rows (); \
00754     octave_idx_type m2_nc = m2.cols (); \
00755     \
00756     if (m1_nr == 1 && m1_nc == 1) \
00757       { \
00758     if (C1 (m1.elem(0,0)) OP C2 (Z2)) \
00759           { \
00760             r = SparseBoolMatrix (m2_nr, m2_nc, true); \
00761             for (octave_idx_type j = 0; j < m2_nc; j++) \
00762               for (octave_idx_type i = m2.cidx(j); i < m2.cidx(j+1); i++) \
00763                 if (! (C1 (m1.elem (0,0)) OP C2 (m2.data(i)))) \
00764                   r.data (m2.ridx (i) + j * m2_nr) = false; \
00765             r.maybe_compress (true); \
00766           } \
00767         else \
00768           { \
00769             r = SparseBoolMatrix (m2_nr, m2_nc, m2.nnz ()); \
00770             r.cidx (0) = static_cast<octave_idx_type> (0); \
00771             octave_idx_type nel = 0; \
00772             for (octave_idx_type j = 0; j < m2_nc; j++) \
00773               { \
00774                 for (octave_idx_type i = m2.cidx(j); i < m2.cidx(j+1); i++) \
00775                   if (C1 (m1.elem (0,0)) OP C2 (m2.data(i))) \
00776                     { \
00777                       r.ridx (nel) = m2.ridx (i); \
00778                       r.data (nel++) = true; \
00779                     } \
00780                 r.cidx (j + 1) = nel; \
00781               } \
00782             r.maybe_compress (false); \
00783           } \
00784       } \
00785     else if (m2_nr == 1 && m2_nc == 1) \
00786       { \
00787         if (C1 (Z1) OP C2 (m2.elem (0,0))) \
00788           { \
00789             r = SparseBoolMatrix (m1_nr, m1_nc, true); \
00790             for (octave_idx_type j = 0; j < m1_nc; j++) \
00791               for (octave_idx_type i = m1.cidx(j); i < m1.cidx(j+1); i++) \
00792                 if (! (C1 (m1.data (i)) OP C2 (m2.elem(0,0)))) \
00793                   r.data (m1.ridx (i) + j * m1_nr) = false; \
00794             r.maybe_compress (true); \
00795           } \
00796         else \
00797           { \
00798             r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz ()); \
00799             r.cidx (0) = static_cast<octave_idx_type> (0); \
00800             octave_idx_type nel = 0; \
00801             for (octave_idx_type j = 0; j < m1_nc; j++) \
00802               { \
00803                 for (octave_idx_type i = m1.cidx(j); i < m1.cidx(j+1); i++) \
00804                   if (C1 (m1.data (i)) OP C2 (m2.elem(0,0))) \
00805                     { \
00806                       r.ridx (nel) = m1.ridx (i); \
00807                       r.data (nel++) = true; \
00808                     } \
00809                 r.cidx (j + 1) = nel; \
00810               } \
00811             r.maybe_compress (false); \
00812           } \
00813       } \
00814     else if (m1_nr == m2_nr && m1_nc == m2_nc) \
00815       { \
00816         if (m1_nr != 0 || m1_nc != 0) \
00817           { \
00818             if (C1 (Z1) OP C2 (Z2)) \
00819               { \
00820                 r = SparseBoolMatrix (m1_nr, m1_nc, true); \
00821                 for (octave_idx_type j = 0; j < m1_nc; j++) \
00822                   { \
00823                      octave_idx_type i1 = m1.cidx (j); \
00824                      octave_idx_type e1 = m1.cidx (j+1); \
00825                      octave_idx_type i2 = m2.cidx (j); \
00826                      octave_idx_type e2 = m2.cidx (j+1); \
00827                      while (i1 < e1 || i2 < e2) \
00828                        { \
00829                          if (i1 == e1 || (i2 < e2 && m1.ridx(i1) > m2.ridx(i2))) \
00830                            { \
00831                              if (! (C1 (Z1) OP C2 (m2.data (i2)))) \
00832                                r.data (m2.ridx (i2) + j * m1_nr) = false; \
00833                              i2++; \
00834                            } \
00835                          else if (i2 == e2 || m1.ridx(i1) < m2.ridx(i2)) \
00836                            { \
00837                              if (! (C1 (m1.data (i1)) OP C2 (Z2))) \
00838                                r.data (m1.ridx (i1) + j * m1_nr) = false; \
00839                              i1++; \
00840                            } \
00841                          else \
00842                            { \
00843                              if (! (C1 (m1.data (i1)) OP C2 (m2.data (i2)))) \
00844                                r.data (m1.ridx (i1) + j * m1_nr) = false; \
00845                              i1++; \
00846                              i2++; \
00847                            } \
00848                        } \
00849                   } \
00850                 r.maybe_compress (true); \
00851               } \
00852             else \
00853               { \
00854                 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz () + m2.nnz ()); \
00855                 r.cidx (0) = static_cast<octave_idx_type> (0); \
00856                 octave_idx_type nel = 0; \
00857                 for (octave_idx_type j = 0; j < m1_nc; j++) \
00858                   { \
00859                      octave_idx_type i1 = m1.cidx (j); \
00860                      octave_idx_type e1 = m1.cidx (j+1); \
00861                      octave_idx_type i2 = m2.cidx (j); \
00862                      octave_idx_type e2 = m2.cidx (j+1); \
00863                      while (i1 < e1 || i2 < e2) \
00864                        { \
00865                          if (i1 == e1 || (i2 < e2 && m1.ridx(i1) > m2.ridx(i2))) \
00866                            { \
00867                              if (C1 (Z1) OP C2 (m2.data (i2))) \
00868                                { \
00869                                  r.ridx (nel) = m2.ridx (i2); \
00870                                  r.data (nel++) = true; \
00871                                } \
00872                              i2++; \
00873                            } \
00874                          else if (i2 == e2 || m1.ridx(i1) < m2.ridx(i2)) \
00875                            { \
00876                              if (C1 (m1.data (i1)) OP C2 (Z2)) \
00877                                { \
00878                                  r.ridx (nel) = m1.ridx (i1); \
00879                                  r.data (nel++) = true; \
00880                                } \
00881                              i1++; \
00882                            } \
00883                          else \
00884                            { \
00885                              if (C1 (m1.data (i1)) OP C2 (m2.data (i2))) \
00886                                { \
00887                                  r.ridx (nel) = m1.ridx (i1); \
00888                                  r.data (nel++) = true; \
00889                                } \
00890                              i1++; \
00891                              i2++; \
00892                            } \
00893                        } \
00894                      r.cidx (j + 1) = nel; \
00895                   } \
00896                 r.maybe_compress (false); \
00897               } \
00898           } \
00899       }       \
00900     else \
00901       { \
00902         if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \
00903           gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
00904       } \
00905     return r; \
00906   }
00907 
00908 #define SPARSE_SMSM_CMP_OPS(M1, Z1, C1, M2, Z2, C2)  \
00909   SPARSE_SMSM_CMP_OP (mx_el_lt, <,  M1, Z1,   , M2, Z2,   ) \
00910   SPARSE_SMSM_CMP_OP (mx_el_le, <=, M1, Z1,   , M2, Z2,   ) \
00911   SPARSE_SMSM_CMP_OP (mx_el_ge, >=, M1, Z1,   , M2, Z2,   ) \
00912   SPARSE_SMSM_CMP_OP (mx_el_gt, >,  M1, Z1,   , M2, Z2,   ) \
00913   SPARSE_SMSM_CMP_OP (mx_el_eq, ==, M1, Z1,   , M2, Z2,   ) \
00914   SPARSE_SMSM_CMP_OP (mx_el_ne, !=, M1, Z1,   , M2, Z2,   )
00915 
00916 #define SPARSE_SMSM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2)  \
00917   SPARSE_SMSM_CMP_OP (mx_el_eq, ==, M1, Z1,   , M2, Z2,   ) \
00918   SPARSE_SMSM_CMP_OP (mx_el_ne, !=, M1, Z1,   , M2, Z2,   )
00919 
00920 #define SPARSE_SMSM_BOOL_OP_DECLS(M1, M2, API) \
00921   SPARSE_BOOL_OP_DECL (mx_el_and, M1, M2, API); \
00922   SPARSE_BOOL_OP_DECL (mx_el_or,  M1, M2, API);
00923 
00924 // FIXME -- this macro duplicatest the bodies of the template
00925 // functions defined in the SPARSE_SSM_BOOL_OP and SPARSE_SMS_BOOL_OP
00926 // macros.
00927 
00928 #define SPARSE_SMSM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \
00929   SparseBoolMatrix \
00930   F (const M1& m1, const M2& m2) \
00931   { \
00932     SparseBoolMatrix r; \
00933     \
00934     octave_idx_type m1_nr = m1.rows (); \
00935     octave_idx_type m1_nc = m1.cols (); \
00936     \
00937     octave_idx_type m2_nr = m2.rows (); \
00938     octave_idx_type m2_nc = m2.cols (); \
00939     \
00940     if (m1_nr == 1 && m1_nc == 1) \
00941       { \
00942         if (m2_nr > 0 && m2_nc > 0) \
00943           { \
00944             if ((m1.elem(0,0) != LHS_ZERO) OP RHS_ZERO) \
00945               { \
00946                 r = SparseBoolMatrix (m2_nr, m2_nc, true); \
00947                 for (octave_idx_type j = 0; j < m2_nc; j++) \
00948                   for (octave_idx_type i = m2.cidx(j); i < m2.cidx(j+1); i++) \
00949                     if (! ((m1.elem(0,0) != LHS_ZERO) OP (m2.data(i) != RHS_ZERO))) \
00950                       r.data (m2.ridx (i) + j * m2_nr) = false; \
00951                 r.maybe_compress (true); \
00952               } \
00953             else \
00954               { \
00955                 r = SparseBoolMatrix (m2_nr, m2_nc, m2.nnz ()); \
00956                 r.cidx (0) = static_cast<octave_idx_type> (0); \
00957                 octave_idx_type nel = 0; \
00958                 for (octave_idx_type j = 0; j < m2_nc; j++) \
00959                   { \
00960                     for (octave_idx_type i = m2.cidx(j); i < m2.cidx(j+1); i++) \
00961                       if ((m1.elem(0,0) != LHS_ZERO) OP (m2.data(i) != RHS_ZERO)) \
00962                         { \
00963                           r.ridx (nel) = m2.ridx (i); \
00964                           r.data (nel++) = true; \
00965                         } \
00966                     r.cidx (j + 1) = nel; \
00967                   } \
00968                 r.maybe_compress (false); \
00969               } \
00970           } \
00971       } \
00972     else if (m2_nr == 1 && m2_nc == 1) \
00973       { \
00974         if (m1_nr > 0 && m1_nc > 0) \
00975           { \
00976             if (LHS_ZERO OP (m2.elem(0,0) != RHS_ZERO)) \
00977               { \
00978                 r = SparseBoolMatrix (m1_nr, m1_nc, true); \
00979                 for (octave_idx_type j = 0; j < m1_nc; j++) \
00980                   for (octave_idx_type i = m1.cidx(j); i < m1.cidx(j+1); i++) \
00981                     if (! ((m1.data(i) != LHS_ZERO) OP (m2.elem(0,0) != RHS_ZERO))) \
00982                       r.data (m1.ridx (i) + j * m1_nr) = false; \
00983                 r.maybe_compress (true); \
00984               } \
00985             else \
00986               { \
00987                 r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz ()); \
00988                 r.cidx (0) = static_cast<octave_idx_type> (0); \
00989                 octave_idx_type nel = 0; \
00990                 for (octave_idx_type j = 0; j < m1_nc; j++) \
00991                   { \
00992                     for (octave_idx_type i = m1.cidx(j); i < m1.cidx(j+1); i++) \
00993                       if ((m1.data(i) != LHS_ZERO) OP (m2.elem(0,0) != RHS_ZERO)) \
00994                         { \
00995                           r.ridx (nel) = m1.ridx (i); \
00996                           r.data (nel++) = true; \
00997                         } \
00998                     r.cidx (j + 1) = nel; \
00999                   } \
01000                 r.maybe_compress (false); \
01001               } \
01002           } \
01003       } \
01004     else if (m1_nr == m2_nr && m1_nc == m2_nc) \
01005       { \
01006         if (m1_nr != 0 || m1_nc != 0) \
01007           { \
01008             r = SparseBoolMatrix (m1_nr, m1_nc, m1.nnz () + m2.nnz ()); \
01009             r.cidx (0) = static_cast<octave_idx_type> (0); \
01010             octave_idx_type nel = 0; \
01011             for (octave_idx_type j = 0; j < m1_nc; j++) \
01012               { \
01013                 octave_idx_type i1 = m1.cidx (j); \
01014                 octave_idx_type e1 = m1.cidx (j+1); \
01015                 octave_idx_type i2 = m2.cidx (j); \
01016                 octave_idx_type e2 = m2.cidx (j+1); \
01017                 while (i1 < e1 || i2 < e2) \
01018                   { \
01019                     if (i1 == e1 || (i2 < e2 && m1.ridx(i1) > m2.ridx(i2))) \
01020                       { \
01021                         if (LHS_ZERO OP m2.data (i2) != RHS_ZERO) \
01022                           { \
01023                             r.ridx (nel) = m2.ridx (i2); \
01024                             r.data (nel++) = true; \
01025                           } \
01026                         i2++; \
01027                       } \
01028                     else if (i2 == e2 || m1.ridx(i1) < m2.ridx(i2)) \
01029                       { \
01030                         if (m1.data (i1) != LHS_ZERO OP RHS_ZERO) \
01031                           { \
01032                             r.ridx (nel) = m1.ridx (i1); \
01033                             r.data (nel++) = true; \
01034                           } \
01035                         i1++; \
01036                       } \
01037                     else \
01038                       { \
01039                         if (m1.data (i1) != LHS_ZERO OP m2.data(i2) != RHS_ZERO) \
01040                           { \
01041                             r.ridx (nel) = m1.ridx (i1); \
01042                             r.data (nel++) = true; \
01043                           } \
01044                         i1++; \
01045                         i2++; \
01046                       } \
01047                   } \
01048                 r.cidx (j + 1) = nel; \
01049               } \
01050             r.maybe_compress (false); \
01051           } \
01052       }       \
01053     else \
01054       { \
01055         if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \
01056           gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
01057       } \
01058     return r; \
01059   }
01060 
01061 #define SPARSE_SMSM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \
01062   SPARSE_SMSM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \
01063   SPARSE_SMSM_BOOL_OP (mx_el_or,  ||, M1, M2, LHS_ZERO, RHS_ZERO) \
01064 
01065 #define SPARSE_SMSM_BOOL_OPS(M1, M2, ZERO) \
01066   SPARSE_SMSM_BOOL_OPS2(M1, M2, ZERO, ZERO)
01067 
01068 #define SPARSE_SMSM_OP_DECLS(R1, R2, M1, M2, API) \
01069   SPARSE_SMSM_BIN_OP_DECLS (R1, R2, M1, M2, API) \
01070   SPARSE_SMSM_CMP_OP_DECLS (M1, M2, API) \
01071   SPARSE_SMSM_BOOL_OP_DECLS (M1, M2, API)
01072 
01073 // matrix by matrix operations.
01074 
01075 #define SPARSE_MSM_BIN_OP_DECLS(R1, R2, M1, M2, API)    \
01076   SPARSE_BIN_OP_DECL (R1, operator +, M1, M2, API); \
01077   SPARSE_BIN_OP_DECL (R1, operator -, M1, M2, API); \
01078   SPARSE_BIN_OP_DECL (R2, product,    M1, M2, API); \
01079   SPARSE_BIN_OP_DECL (R2, quotient,   M1, M2, API);
01080 
01081 #define SPARSE_MSM_BIN_OP_1(R, F, OP, M1, M2)   \
01082   R \
01083   F (const M1& m1, const M2& m2) \
01084   { \
01085     R r; \
01086  \
01087     octave_idx_type m1_nr = m1.rows (); \
01088     octave_idx_type m1_nc = m1.cols (); \
01089  \
01090     octave_idx_type m2_nr = m2.rows (); \
01091     octave_idx_type m2_nc = m2.cols (); \
01092  \
01093     if (m2_nr == 1 && m2_nc == 1) \
01094       r = R (m1 OP m2.elem(0,0)); \
01095     else if (m1_nr != m2_nr || m1_nc != m2_nc) \
01096       gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
01097     else \
01098       { \
01099         r = R (F (m1, m2.matrix_value ())); \
01100       } \
01101     return r; \
01102   }
01103 
01104 #define SPARSE_MSM_BIN_OP_2(R, F, OP, M1, M2) \
01105   R \
01106   F (const M1& m1, const M2& m2) \
01107   { \
01108     R r; \
01109  \
01110     octave_idx_type m1_nr = m1.rows (); \
01111     octave_idx_type m1_nc = m1.cols (); \
01112  \
01113     octave_idx_type m2_nr = m2.rows (); \
01114     octave_idx_type m2_nc = m2.cols (); \
01115  \
01116     if (m2_nr == 1 && m2_nc == 1) \
01117       r = R (m1 OP m2.elem(0,0)); \
01118     else if (m1_nr != m2_nr || m1_nc != m2_nc) \
01119       gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
01120     else \
01121       { \
01122         if (do_mx_check (m1, mx_inline_all_finite<M1::element_type>)) \
01123           { \
01124             /* Sparsity pattern is preserved. */ \
01125             octave_idx_type m2_nz = m2.nnz (); \
01126             r = R (m2_nr, m2_nc, m2_nz); \
01127             for (octave_idx_type j = 0, k = 0; j < m2_nc; j++) \
01128               { \
01129                 octave_quit (); \
01130                 for (octave_idx_type i = m2.cidx(j); i < m2.cidx(j+1); i++) \
01131                   { \
01132                     octave_idx_type mri = m2.ridx(i); \
01133                     R::element_type x = m1(mri, j) OP m2.data(i); \
01134                     if (x != 0.0) \
01135                       { \
01136                         r.xdata(k) = x; \
01137                         r.xridx(k) = m2.ridx(i); \
01138                         k++; \
01139                       } \
01140                   } \
01141                 r.xcidx(j+1) = k; \
01142               } \
01143             r.maybe_compress (false); \
01144             return r; \
01145           } \
01146         else \
01147           r = R (F (m1, m2.matrix_value ())); \
01148       } \
01149  \
01150     return r; \
01151   }
01152 
01153 // FIXME Pass a specific ZERO value
01154 #define SPARSE_MSM_BIN_OPS(R1, R2, M1, M2) \
01155   SPARSE_MSM_BIN_OP_1 (R1, operator +,  +, M1, M2) \
01156   SPARSE_MSM_BIN_OP_1 (R1, operator -,  -, M1, M2) \
01157   SPARSE_MSM_BIN_OP_2 (R2, product,     *, M1, M2) \
01158   SPARSE_MSM_BIN_OP_1 (R2, quotient,    /, M1, M2)
01159 
01160 #define SPARSE_MSM_CMP_OP_DECLS(M1, M2, API) \
01161   SPARSE_CMP_OP_DECL (mx_el_lt, M1, M2, API); \
01162   SPARSE_CMP_OP_DECL (mx_el_le, M1, M2, API); \
01163   SPARSE_CMP_OP_DECL (mx_el_ge, M1, M2, API); \
01164   SPARSE_CMP_OP_DECL (mx_el_gt, M1, M2, API); \
01165   SPARSE_CMP_OP_DECL (mx_el_eq, M1, M2, API); \
01166   SPARSE_CMP_OP_DECL (mx_el_ne, M1, M2, API);
01167 
01168 #define SPARSE_MSM_EQNE_OP_DECLS(M1, M2, API) \
01169   SPARSE_CMP_OP_DECL (mx_el_eq, M1, M2, API); \
01170   SPARSE_CMP_OP_DECL (mx_el_ne, M1, M2, API);
01171 
01172 #define SPARSE_MSM_CMP_OP(F, OP, M1, C1, M2, C2)        \
01173   SparseBoolMatrix \
01174   F (const M1& m1, const M2& m2) \
01175   { \
01176     SparseBoolMatrix r; \
01177     \
01178     octave_idx_type m1_nr = m1.rows (); \
01179     octave_idx_type m1_nc = m1.cols (); \
01180     \
01181     octave_idx_type m2_nr = m2.rows (); \
01182     octave_idx_type m2_nc = m2.cols (); \
01183     \
01184     if (m2_nr == 1 && m2_nc == 1) \
01185       r = SparseBoolMatrix (F (m1, m2.elem(0,0))); \
01186     else if (m1_nr == m2_nr && m1_nc == m2_nc) \
01187       { \
01188         if (m1_nr != 0 || m1_nc != 0) \
01189           { \
01190             /* Count num of non-zero elements */ \
01191             octave_idx_type nel = 0; \
01192             for (octave_idx_type j = 0; j < m1_nc; j++) \
01193               for (octave_idx_type i = 0; i < m1_nr; i++) \
01194                 if (C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j))) \
01195                   nel++; \
01196             \
01197             r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
01198             \
01199             octave_idx_type ii = 0; \
01200             r.cidx (0) = 0; \
01201             for (octave_idx_type j = 0; j < m1_nc; j++) \
01202               { \
01203                 for (octave_idx_type i = 0; i < m1_nr; i++) \
01204                   { \
01205                     bool el = C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j)); \
01206                     if (el) \
01207                       { \
01208                         r.data(ii) = el; \
01209                         r.ridx(ii++) = i; \
01210                       } \
01211                   } \
01212                 r.cidx(j+1) = ii; \
01213               } \
01214           } \
01215       }       \
01216     else \
01217       { \
01218         if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \
01219           gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
01220       } \
01221     return r; \
01222   }
01223 
01224 #define SPARSE_MSM_CMP_OPS(M1, Z1, C1, M2, Z2, C2)  \
01225   SPARSE_MSM_CMP_OP (mx_el_lt, <,  M1,   , M2,   ) \
01226   SPARSE_MSM_CMP_OP (mx_el_le, <=, M1,   , M2,   ) \
01227   SPARSE_MSM_CMP_OP (mx_el_ge, >=, M1,   , M2,   ) \
01228   SPARSE_MSM_CMP_OP (mx_el_gt, >,  M1,   , M2,   ) \
01229   SPARSE_MSM_CMP_OP (mx_el_eq, ==, M1,   , M2,   ) \
01230   SPARSE_MSM_CMP_OP (mx_el_ne, !=, M1,   , M2,   )
01231 
01232 #define SPARSE_MSM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2)  \
01233   SPARSE_MSM_CMP_OP (mx_el_eq, ==, M1,   , M2,   ) \
01234   SPARSE_MSM_CMP_OP (mx_el_ne, !=, M1,   , M2,   )
01235 
01236 #define SPARSE_MSM_BOOL_OP_DECLS(M1, M2, API) \
01237   SPARSE_BOOL_OP_DECL (mx_el_and, M1, M2, API); \
01238   SPARSE_BOOL_OP_DECL (mx_el_or,  M1, M2, API);
01239 
01240 #define SPARSE_MSM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \
01241   SparseBoolMatrix \
01242   F (const M1& m1, const M2& m2) \
01243   { \
01244     SparseBoolMatrix r; \
01245     \
01246     octave_idx_type m1_nr = m1.rows (); \
01247     octave_idx_type m1_nc = m1.cols (); \
01248     \
01249     octave_idx_type m2_nr = m2.rows (); \
01250     octave_idx_type m2_nc = m2.cols (); \
01251     \
01252     if (m2_nr == 1 && m2_nc == 1) \
01253       r = SparseBoolMatrix  (F (m1, m2.elem(0,0))); \
01254     else if (m1_nr == m2_nr && m1_nc == m2_nc) \
01255       { \
01256         if (m1_nr != 0 || m1_nc != 0) \
01257           { \
01258             /* Count num of non-zero elements */ \
01259             octave_idx_type nel = 0; \
01260             for (octave_idx_type j = 0; j < m1_nc; j++) \
01261               for (octave_idx_type i = 0; i < m1_nr; i++) \
01262                 if ((m1.elem(i, j) != LHS_ZERO) \
01263                     OP (m2.elem(i, j) != RHS_ZERO)) \
01264                   nel++; \
01265             \
01266             r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
01267             \
01268             octave_idx_type ii = 0; \
01269             r.cidx (0) = 0; \
01270             for (octave_idx_type j = 0; j < m1_nc; j++) \
01271               { \
01272                 for (octave_idx_type i = 0; i < m1_nr; i++) \
01273                   { \
01274                     bool el = (m1.elem(i, j) != LHS_ZERO) \
01275                       OP (m2.elem(i, j) != RHS_ZERO);     \
01276                     if (el) \
01277                       { \
01278                         r.data(ii) = el; \
01279                         r.ridx(ii++) = i; \
01280                       } \
01281                   } \
01282                 r.cidx(j+1) = ii; \
01283               } \
01284           } \
01285       }       \
01286     else \
01287       { \
01288         if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \
01289           gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
01290       } \
01291     return r; \
01292   }
01293 
01294 #define SPARSE_MSM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \
01295   SPARSE_MSM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \
01296   SPARSE_MSM_BOOL_OP (mx_el_or,  ||, M1, M2, LHS_ZERO, RHS_ZERO) \
01297 
01298 #define SPARSE_MSM_BOOL_OPS(M1, M2, ZERO) \
01299   SPARSE_MSM_BOOL_OPS2(M1, M2, ZERO, ZERO)
01300 
01301 #define SPARSE_MSM_OP_DECLS(R1, R2, M1, M2, API) \
01302   SPARSE_MSM_BIN_OP_DECLS (R1, R2, M1, M2, API) \
01303   SPARSE_MSM_CMP_OP_DECLS (M1, M2, API) \
01304   SPARSE_MSM_BOOL_OP_DECLS (M1, M2, API)
01305 
01306 // matrix by matrix operations.
01307 
01308 #define SPARSE_SMM_BIN_OP_DECLS(R1, R2, M1, M2, API)    \
01309   SPARSE_BIN_OP_DECL (R1, operator +, M1, M2, API); \
01310   SPARSE_BIN_OP_DECL (R1, operator -, M1, M2, API); \
01311   SPARSE_BIN_OP_DECL (R2, product,    M1, M2, API); \
01312   SPARSE_BIN_OP_DECL (R2, quotient,   M1, M2, API);
01313 
01314 #define SPARSE_SMM_BIN_OP_1(R, F, OP, M1, M2)   \
01315   R \
01316   F (const M1& m1, const M2& m2) \
01317   { \
01318     R r; \
01319  \
01320     octave_idx_type m1_nr = m1.rows (); \
01321     octave_idx_type m1_nc = m1.cols (); \
01322  \
01323     octave_idx_type m2_nr = m2.rows (); \
01324     octave_idx_type m2_nc = m2.cols (); \
01325  \
01326     if (m1_nr == 1 && m1_nc == 1) \
01327       r = R (m1.elem(0,0) OP m2); \
01328     else if (m1_nr != m2_nr || m1_nc != m2_nc) \
01329       gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
01330     else \
01331       { \
01332         r = R (m1.matrix_value () OP m2); \
01333       } \
01334     return r; \
01335   }
01336 
01337 // sm .* m preserves sparsity if m contains no Infs nor Nans.
01338 #define SPARSE_SMM_BIN_OP_2_CHECK_product(ET) \
01339   do_mx_check (m2, mx_inline_all_finite<ET>)
01340 
01341 // sm ./ m preserves sparsity if m contains no NaNs or zeros.
01342 #define SPARSE_SMM_BIN_OP_2_CHECK_quotient(ET) \
01343   ! do_mx_check (m2, mx_inline_any_nan<ET>) && m2.nnz () == m2.numel ()
01344 
01345 #define SPARSE_SMM_BIN_OP_2(R, F, OP, M1, M2) \
01346   R \
01347   F (const M1& m1, const M2& m2) \
01348   { \
01349     R r; \
01350  \
01351     octave_idx_type m1_nr = m1.rows (); \
01352     octave_idx_type m1_nc = m1.cols (); \
01353  \
01354     octave_idx_type m2_nr = m2.rows (); \
01355     octave_idx_type m2_nc = m2.cols (); \
01356  \
01357     if (m1_nr == 1 && m1_nc == 1) \
01358       r = R (m1.elem(0,0) OP m2); \
01359     else if (m1_nr != m2_nr || m1_nc != m2_nc) \
01360       gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
01361     else \
01362       { \
01363         if (SPARSE_SMM_BIN_OP_2_CHECK_ ## F(M2::element_type)) \
01364           { \
01365             /* Sparsity pattern is preserved. */ \
01366             octave_idx_type m1_nz = m1.nnz (); \
01367             r = R (m1_nr, m1_nc, m1_nz); \
01368             for (octave_idx_type j = 0, k = 0; j < m1_nc; j++) \
01369               { \
01370                 octave_quit (); \
01371                 for (octave_idx_type i = m1.cidx(j); i < m1.cidx(j+1); i++) \
01372                   { \
01373                     octave_idx_type mri = m1.ridx(i); \
01374                     R::element_type x = m1.data(i) OP m2(mri, j); \
01375                     if (x != 0.0) \
01376                       { \
01377                         r.xdata(k) = x; \
01378                         r.xridx(k) = m1.ridx(i); \
01379                         k++; \
01380                       } \
01381                   } \
01382                 r.xcidx(j+1) = k; \
01383               } \
01384             r.maybe_compress (false); \
01385             return r; \
01386           } \
01387         else \
01388           r = R (F (m1.matrix_value (), m2)); \
01389       } \
01390  \
01391     return r; \
01392   }
01393 
01394 #define SPARSE_SMM_BIN_OPS(R1, R2, M1, M2) \
01395   SPARSE_SMM_BIN_OP_1 (R1, operator +,  +, M1, M2) \
01396   SPARSE_SMM_BIN_OP_1 (R1, operator -,  -, M1, M2) \
01397   SPARSE_SMM_BIN_OP_2 (R2, product,     *, M1, M2) \
01398   SPARSE_SMM_BIN_OP_2 (R2, quotient,    /, M1, M2)
01399 
01400 #define SPARSE_SMM_CMP_OP_DECLS(M1, M2, API) \
01401   SPARSE_CMP_OP_DECL (mx_el_lt, M1, M2, API); \
01402   SPARSE_CMP_OP_DECL (mx_el_le, M1, M2, API); \
01403   SPARSE_CMP_OP_DECL (mx_el_ge, M1, M2, API); \
01404   SPARSE_CMP_OP_DECL (mx_el_gt, M1, M2, API); \
01405   SPARSE_CMP_OP_DECL (mx_el_eq, M1, M2, API); \
01406   SPARSE_CMP_OP_DECL (mx_el_ne, M1, M2, API);
01407 
01408 #define SPARSE_SMM_EQNE_OP_DECLS(M1, M2, API) \
01409   SPARSE_CMP_OP_DECL (mx_el_eq, M1, M2, API); \
01410   SPARSE_CMP_OP_DECL (mx_el_ne, M1, M2, API);
01411 
01412 #define SPARSE_SMM_CMP_OP(F, OP, M1, C1, M2, C2)        \
01413   SparseBoolMatrix \
01414   F (const M1& m1, const M2& m2) \
01415   { \
01416     SparseBoolMatrix r; \
01417     \
01418     octave_idx_type m1_nr = m1.rows (); \
01419     octave_idx_type m1_nc = m1.cols (); \
01420     \
01421     octave_idx_type m2_nr = m2.rows (); \
01422     octave_idx_type m2_nc = m2.cols (); \
01423     \
01424     if (m1_nr == 1 && m1_nc == 1) \
01425       r = SparseBoolMatrix (F (m1.elem(0,0), m2)); \
01426     else if (m1_nr == m2_nr && m1_nc == m2_nc) \
01427       { \
01428         if (m1_nr != 0 || m1_nc != 0) \
01429           { \
01430             /* Count num of non-zero elements */ \
01431             octave_idx_type nel = 0; \
01432             for (octave_idx_type j = 0; j < m1_nc; j++) \
01433               for (octave_idx_type i = 0; i < m1_nr; i++) \
01434                 if (C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j))) \
01435                   nel++; \
01436             \
01437             r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
01438             \
01439             octave_idx_type ii = 0; \
01440             r.cidx (0) = 0; \
01441             for (octave_idx_type j = 0; j < m1_nc; j++) \
01442               { \
01443                 for (octave_idx_type i = 0; i < m1_nr; i++) \
01444                   { \
01445                     bool el = C1 (m1.elem(i, j)) OP C2 (m2.elem(i, j)); \
01446                     if (el) \
01447                       { \
01448                         r.data(ii) = el; \
01449                         r.ridx(ii++) = i; \
01450                       } \
01451                   } \
01452                 r.cidx(j+1) = ii; \
01453               } \
01454           } \
01455       }       \
01456     else \
01457       { \
01458         if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \
01459           gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
01460       } \
01461     return r; \
01462   }
01463 
01464 #define SPARSE_SMM_CMP_OPS(M1, Z1, C1, M2, Z2, C2)  \
01465   SPARSE_SMM_CMP_OP (mx_el_lt, <,  M1,   , M2,   ) \
01466   SPARSE_SMM_CMP_OP (mx_el_le, <=, M1,   , M2,   ) \
01467   SPARSE_SMM_CMP_OP (mx_el_ge, >=, M1,   , M2,   ) \
01468   SPARSE_SMM_CMP_OP (mx_el_gt, >,  M1,   , M2,   ) \
01469   SPARSE_SMM_CMP_OP (mx_el_eq, ==, M1,   , M2,   ) \
01470   SPARSE_SMM_CMP_OP (mx_el_ne, !=, M1,   , M2,   )
01471 
01472 #define SPARSE_SMM_EQNE_OPS(M1, Z1, C1, M2, Z2, C2)  \
01473   SPARSE_SMM_CMP_OP (mx_el_eq, ==, M1,   , M2,   ) \
01474   SPARSE_SMM_CMP_OP (mx_el_ne, !=, M1,   , M2,   )
01475 
01476 #define SPARSE_SMM_BOOL_OP_DECLS(M1, M2, API) \
01477   SPARSE_BOOL_OP_DECL (mx_el_and, M1, M2, API); \
01478   SPARSE_BOOL_OP_DECL (mx_el_or,  M1, M2, API);
01479 
01480 #define SPARSE_SMM_BOOL_OP(F, OP, M1, M2, LHS_ZERO, RHS_ZERO) \
01481   SparseBoolMatrix \
01482   F (const M1& m1, const M2& m2) \
01483   { \
01484     SparseBoolMatrix r; \
01485     \
01486     octave_idx_type m1_nr = m1.rows (); \
01487     octave_idx_type m1_nc = m1.cols (); \
01488     \
01489     octave_idx_type m2_nr = m2.rows (); \
01490     octave_idx_type m2_nc = m2.cols (); \
01491     \
01492     if (m1_nr == 1 && m1_nc == 1) \
01493       r = SparseBoolMatrix (F (m1.elem(0,0), m2)); \
01494     else if (m1_nr == m2_nr && m1_nc == m2_nc) \
01495       { \
01496         if (m1_nr != 0 || m1_nc != 0) \
01497           { \
01498             /* Count num of non-zero elements */ \
01499             octave_idx_type nel = 0; \
01500             for (octave_idx_type j = 0; j < m1_nc; j++) \
01501               for (octave_idx_type i = 0; i < m1_nr; i++) \
01502                 if ((m1.elem(i, j) != LHS_ZERO) \
01503                     OP (m2.elem(i, j) != RHS_ZERO)) \
01504                   nel++; \
01505             \
01506             r = SparseBoolMatrix (m1_nr, m1_nc, nel); \
01507             \
01508             octave_idx_type ii = 0; \
01509             r.cidx (0) = 0; \
01510             for (octave_idx_type j = 0; j < m1_nc; j++) \
01511               { \
01512                 for (octave_idx_type i = 0; i < m1_nr; i++) \
01513                   { \
01514                     bool el = (m1.elem(i, j) != LHS_ZERO) \
01515                       OP (m2.elem(i, j) != RHS_ZERO);     \
01516                     if (el) \
01517                       { \
01518                         r.data(ii) = el; \
01519                         r.ridx(ii++) = i; \
01520                       } \
01521                   } \
01522                 r.cidx(j+1) = ii; \
01523               } \
01524           } \
01525       }       \
01526     else \
01527       { \
01528         if ((m1_nr != 0 || m1_nc != 0) && (m2_nr != 0 || m2_nc != 0)) \
01529           gripe_nonconformant (#F, m1_nr, m1_nc, m2_nr, m2_nc); \
01530       } \
01531     return r; \
01532   }
01533 
01534 #define SPARSE_SMM_BOOL_OPS2(M1, M2, LHS_ZERO, RHS_ZERO) \
01535   SPARSE_SMM_BOOL_OP (mx_el_and, &&, M1, M2, LHS_ZERO, RHS_ZERO) \
01536   SPARSE_SMM_BOOL_OP (mx_el_or,  ||, M1, M2, LHS_ZERO, RHS_ZERO) \
01537 
01538 #define SPARSE_SMM_BOOL_OPS(M1, M2, ZERO) \
01539   SPARSE_SMM_BOOL_OPS2(M1, M2, ZERO, ZERO)
01540 
01541 #define SPARSE_SMM_OP_DECLS(R1, R2, M1, M2, API) \
01542   SPARSE_SMM_BIN_OP_DECLS (R1, R2, M1, M2, API) \
01543   SPARSE_SMM_CMP_OP_DECLS (M1, M2, API) \
01544   SPARSE_SMM_BOOL_OP_DECLS (M1, M2, API)
01545 
01546 // Avoid some code duplication.  Maybe we should use templates.
01547 
01548 #define SPARSE_CUMSUM(RET_TYPE, ELT_TYPE, FCN)  \
01549  \
01550   octave_idx_type nr = rows (); \
01551   octave_idx_type nc = cols (); \
01552  \
01553   RET_TYPE retval; \
01554  \
01555   if (nr > 0 && nc > 0) \
01556     { \
01557       if ((nr == 1 && dim == -1) || dim == 1) \
01558         /* Ugly!! Is there a better way? */ \
01559         retval = transpose (). FCN (0) .transpose (); \
01560       else \
01561         { \
01562           octave_idx_type nel = 0; \
01563           for (octave_idx_type i = 0; i < nc; i++) \
01564             { \
01565               ELT_TYPE t = ELT_TYPE (); \
01566               for (octave_idx_type j = cidx (i); j < cidx (i+1); j++)   \
01567                 { \
01568                   t += data(j); \
01569                   if (t != ELT_TYPE ()) \
01570                     { \
01571                       if (j == cidx(i+1) - 1) \
01572                         nel += nr - ridx(j);  \
01573                       else \
01574                         nel += ridx(j+1) - ridx(j); \
01575                     } \
01576                 } \
01577             } \
01578           retval = RET_TYPE (nr, nc, nel); \
01579           retval.cidx(0) = 0; \
01580           octave_idx_type ii = 0; \
01581           for (octave_idx_type i = 0; i < nc; i++) \
01582             { \
01583               ELT_TYPE t = ELT_TYPE (); \
01584               for (octave_idx_type j = cidx (i); j < cidx (i+1); j++)   \
01585                 { \
01586                   t += data(j); \
01587                   if (t != ELT_TYPE ()) \
01588                     { \
01589                       if (j == cidx(i+1) - 1) \
01590                         { \
01591                           for (octave_idx_type k = ridx(j); k < nr; k++) \
01592                             { \
01593                                retval.data (ii) = t; \
01594                                retval.ridx (ii++) = k; \
01595                             } \
01596                         } \
01597                       else \
01598                         { \
01599                           for (octave_idx_type k = ridx(j); k < ridx(j+1); k++) \
01600                             { \
01601                                retval.data (ii) = t; \
01602                                retval.ridx (ii++) = k; \
01603                             } \
01604                         } \
01605                     } \
01606                 } \
01607               retval.cidx(i+1) = ii; \
01608             } \
01609         } \
01610     } \
01611   else \
01612     retval = RET_TYPE (nr,nc); \
01613  \
01614   return retval
01615 
01616 
01617 #define SPARSE_CUMPROD(RET_TYPE, ELT_TYPE, FCN) \
01618  \
01619   octave_idx_type nr = rows (); \
01620   octave_idx_type nc = cols (); \
01621  \
01622   RET_TYPE retval; \
01623  \
01624   if (nr > 0 && nc > 0) \
01625     { \
01626       if ((nr == 1 && dim == -1) || dim == 1) \
01627         /* Ugly!! Is there a better way? */ \
01628         retval = transpose (). FCN (0) .transpose (); \
01629       else \
01630         { \
01631           octave_idx_type nel = 0; \
01632           for (octave_idx_type i = 0; i < nc; i++) \
01633             { \
01634               octave_idx_type jj = 0; \
01635               for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
01636                 { \
01637                   if (jj == ridx(j)) \
01638                     { \
01639                       nel++; \
01640                       jj++; \
01641                     } \
01642                   else \
01643                     break; \
01644                 } \
01645             } \
01646           retval = RET_TYPE (nr, nc, nel); \
01647           retval.cidx(0) = 0; \
01648           octave_idx_type ii = 0; \
01649           for (octave_idx_type i = 0; i < nc; i++) \
01650             { \
01651               ELT_TYPE t = ELT_TYPE (1.); \
01652               octave_idx_type jj = 0; \
01653               for (octave_idx_type j = cidx (i); j < cidx (i+1); j++) \
01654                 { \
01655                   if (jj == ridx(j)) \
01656                     { \
01657                       t *= data(j); \
01658                       retval.data(ii) = t; \
01659                       retval.ridx(ii++) = jj++; \
01660                     } \
01661                   else \
01662                     break; \
01663                 } \
01664               retval.cidx(i+1) = ii; \
01665             } \
01666         } \
01667     } \
01668   else \
01669     retval = RET_TYPE (nr,nc); \
01670  \
01671   return retval
01672 
01673 #define SPARSE_BASE_REDUCTION_OP(RET_TYPE, EL_TYPE, ROW_EXPR, COL_EXPR, \
01674                                  INIT_VAL, MT_RESULT) \
01675  \
01676   octave_idx_type nr = rows (); \
01677   octave_idx_type nc = cols (); \
01678  \
01679   RET_TYPE retval; \
01680  \
01681   if (nr > 0 && nc > 0) \
01682     { \
01683       if ((nr == 1 && dim == -1) || dim == 1) \
01684         { \
01685           /* Define j here to allow fancy definition for prod method */ \
01686           octave_idx_type j = 0; \
01687           OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nr); \
01688           \
01689           for (octave_idx_type i = 0; i < nr; i++) \
01690             tmp[i] = INIT_VAL; \
01691           for (j = 0; j < nc; j++) \
01692             { \
01693               for (octave_idx_type i = cidx(j); i < cidx(j + 1); i++) \
01694                 { \
01695                   ROW_EXPR; \
01696                 } \
01697             } \
01698           octave_idx_type nel = 0; \
01699           for (octave_idx_type i = 0; i < nr; i++) \
01700             if (tmp[i] != EL_TYPE ())  \
01701               nel++ ; \
01702           retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), nel); \
01703           retval.cidx(0) = 0; \
01704           retval.cidx(1) = nel; \
01705           nel = 0; \
01706           for (octave_idx_type i = 0; i < nr; i++) \
01707             if (tmp[i] != EL_TYPE ())  \
01708               { \
01709                 retval.data(nel) = tmp[i]; \
01710                 retval.ridx(nel++) = i; \
01711               } \
01712         } \
01713       else \
01714         { \
01715           OCTAVE_LOCAL_BUFFER (EL_TYPE, tmp, nc); \
01716           \
01717           for (octave_idx_type j = 0; j < nc; j++) \
01718             { \
01719               tmp[j] = INIT_VAL; \
01720               for (octave_idx_type i = cidx(j); i < cidx(j + 1); i++) \
01721                 { \
01722                   COL_EXPR; \
01723                 } \
01724             } \
01725           octave_idx_type nel = 0; \
01726           for (octave_idx_type i = 0; i < nc; i++) \
01727             if (tmp[i] != EL_TYPE ())  \
01728               nel++ ; \
01729           retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, nel); \
01730           retval.cidx(0) = 0; \
01731           nel = 0; \
01732           for (octave_idx_type i = 0; i < nc; i++) \
01733             if (tmp[i] != EL_TYPE ())  \
01734               { \
01735                 retval.data(nel) = tmp[i]; \
01736                 retval.ridx(nel++) = 0; \
01737                 retval.cidx(i+1) = retval.cidx(i) + 1; \
01738               } \
01739             else \
01740               retval.cidx(i+1) = retval.cidx(i); \
01741         } \
01742     } \
01743   else if (nc == 0 && (nr == 0 || (nr == 1 && dim == -1))) \
01744     { \
01745       if (MT_RESULT) \
01746         { \
01747           retval = RET_TYPE (static_cast<octave_idx_type> (1), \
01748                              static_cast<octave_idx_type> (1), \
01749                              static_cast<octave_idx_type> (1)); \
01750           retval.cidx(0) = 0; \
01751           retval.cidx(1) = 1; \
01752           retval.ridx(0) = 0; \
01753           retval.data(0) = MT_RESULT; \
01754         } \
01755       else \
01756           retval = RET_TYPE (static_cast<octave_idx_type> (1), \
01757                              static_cast<octave_idx_type> (1), \
01758                              static_cast<octave_idx_type> (0)); \
01759     } \
01760   else if (nr == 0 && (dim == 0 || dim == -1)) \
01761     { \
01762       if (MT_RESULT) \
01763         { \
01764           retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, nc); \
01765           retval.cidx (0) = 0; \
01766           for (octave_idx_type i = 0; i < nc ; i++) \
01767             { \
01768               retval.ridx (i) = 0; \
01769               retval.cidx (i+1) = i; \
01770               retval.data (i) = MT_RESULT; \
01771             } \
01772         } \
01773       else \
01774         retval = RET_TYPE (static_cast<octave_idx_type> (1), nc, \
01775                            static_cast<octave_idx_type> (0)); \
01776     } \
01777   else if (nc == 0 && dim == 1) \
01778     { \
01779       if (MT_RESULT) \
01780         { \
01781           retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), nr); \
01782           retval.cidx(0) = 0; \
01783           retval.cidx(1) = nr; \
01784           for (octave_idx_type i = 0; i < nr; i++) \
01785             { \
01786               retval.ridx(i) = i; \
01787               retval.data(i) = MT_RESULT; \
01788             } \
01789         } \
01790       else \
01791         retval = RET_TYPE (nr, static_cast<octave_idx_type> (1), \
01792                            static_cast<octave_idx_type> (0)); \
01793     } \
01794   else \
01795     retval.resize (nr > 0, nc > 0); \
01796  \
01797   return retval
01798 
01799 #define SPARSE_REDUCTION_OP_ROW_EXPR(OP) \
01800   tmp[ridx(i)] OP data (i)
01801 
01802 #define SPARSE_REDUCTION_OP_COL_EXPR(OP) \
01803   tmp[j] OP data (i)
01804 
01805 #define SPARSE_REDUCTION_OP(RET_TYPE, EL_TYPE, OP, INIT_VAL, MT_RESULT) \
01806   SPARSE_BASE_REDUCTION_OP (RET_TYPE, EL_TYPE, \
01807                         SPARSE_REDUCTION_OP_ROW_EXPR (OP), \
01808                         SPARSE_REDUCTION_OP_COL_EXPR (OP), \
01809                         INIT_VAL, MT_RESULT)
01810 
01811 
01812 // Don't break from this loop if the test succeeds because
01813 // we are looping over the rows and not the columns in the inner
01814 // loop.
01815 #define SPARSE_ANY_ALL_OP_ROW_CODE(TEST_OP, TEST_TRUE_VAL) \
01816   if (data (i) TEST_OP 0.0) \
01817     tmp[ridx(i)] = TEST_TRUE_VAL; \
01818 
01819 #define SPARSE_ANY_ALL_OP_COL_CODE(TEST_OP, TEST_TRUE_VAL) \
01820   if (data (i) TEST_OP 0.0) \
01821     { \
01822       tmp[j] = TEST_TRUE_VAL; \
01823       break; \
01824     }
01825 
01826 #define SPARSE_ANY_ALL_OP(DIM, INIT_VAL, MT_RESULT, TEST_OP, TEST_TRUE_VAL) \
01827   SPARSE_BASE_REDUCTION_OP (SparseBoolMatrix, char, \
01828                         SPARSE_ANY_ALL_OP_ROW_CODE (TEST_OP, TEST_TRUE_VAL), \
01829                         SPARSE_ANY_ALL_OP_COL_CODE (TEST_OP, TEST_TRUE_VAL), \
01830                         INIT_VAL, MT_RESULT)
01831 
01832 #define SPARSE_ALL_OP(DIM) \
01833   if ((rows() == 1 && dim == -1) || dim == 1) \
01834     return transpose (). all (0). transpose(); \
01835   else \
01836     { \
01837       SPARSE_ANY_ALL_OP (DIM, (cidx(j+1) - cidx(j) < nr ? false : true), \
01838                          true, ==, false); \
01839     }
01840 
01841 #define SPARSE_ANY_OP(DIM) SPARSE_ANY_ALL_OP (DIM, false, false, !=, true)
01842 
01843 #define SPARSE_SPARSE_MUL( RET_TYPE, RET_EL_TYPE, EL_TYPE ) \
01844   octave_idx_type nr = m.rows (); \
01845   octave_idx_type nc = m.cols (); \
01846   \
01847   octave_idx_type a_nr = a.rows (); \
01848   octave_idx_type a_nc = a.cols (); \
01849   \
01850   if (nr == 1 && nc == 1) \
01851    { \
01852      RET_EL_TYPE s = m.elem(0,0); \
01853      octave_idx_type nz = a.nnz(); \
01854      RET_TYPE r (a_nr, a_nc, nz); \
01855      \
01856      for (octave_idx_type i = 0; i < nz; i++) \
01857        { \
01858          octave_quit (); \
01859          r.data(i) = s * a.data(i); \
01860          r.ridx(i) = a.ridx(i); \
01861        } \
01862      for (octave_idx_type i = 0; i < a_nc + 1; i++) \
01863        { \
01864          octave_quit (); \
01865          r.cidx(i) = a.cidx(i); \
01866        } \
01867      \
01868      r.maybe_compress (true); \
01869      return r; \
01870    } \
01871   else if (a_nr == 1 && a_nc == 1) \
01872    { \
01873      RET_EL_TYPE s = a.elem(0,0); \
01874      octave_idx_type nz = m.nnz(); \
01875      RET_TYPE r (nr, nc, nz); \
01876      \
01877      for (octave_idx_type i = 0; i < nz; i++) \
01878        { \
01879          octave_quit (); \
01880          r.data(i) = m.data(i) * s; \
01881          r.ridx(i) = m.ridx(i); \
01882        } \
01883      for (octave_idx_type i = 0; i < nc + 1; i++) \
01884        { \
01885          octave_quit (); \
01886          r.cidx(i) = m.cidx(i); \
01887        } \
01888      \
01889      r.maybe_compress (true); \
01890      return r; \
01891    } \
01892   else if (nc != a_nr) \
01893     { \
01894       gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc); \
01895       return RET_TYPE (); \
01896     } \
01897   else \
01898     { \
01899       OCTAVE_LOCAL_BUFFER (octave_idx_type, w, nr); \
01900       RET_TYPE retval (nr, a_nc, static_cast<octave_idx_type> (0)); \
01901       for (octave_idx_type i = 0; i < nr; i++) \
01902         w[i] = 0; \
01903       retval.xcidx(0) = 0; \
01904       \
01905       octave_idx_type nel = 0; \
01906       \
01907       for (octave_idx_type i = 0; i < a_nc; i++) \
01908         { \
01909           for (octave_idx_type j = a.cidx(i); j < a.cidx(i+1); j++) \
01910             { \
01911               octave_idx_type  col = a.ridx(j); \
01912               for (octave_idx_type k = m.cidx(col) ; k < m.cidx(col+1); k++) \
01913                 { \
01914                   if (w[m.ridx(k)] < i + 1) \
01915                     { \
01916                       w[m.ridx(k)] = i + 1; \
01917                       nel++; \
01918                     } \
01919                   octave_quit (); \
01920                 } \
01921             } \
01922           retval.xcidx(i+1) = nel; \
01923         } \
01924       \
01925       if (nel == 0) \
01926         return RET_TYPE (nr, a_nc); \
01927       else \
01928         {  \
01929           for (octave_idx_type i = 0; i < nr; i++) \
01930             w[i] = 0; \
01931           \
01932           OCTAVE_LOCAL_BUFFER (RET_EL_TYPE, Xcol, nr); \
01933           \
01934           retval.change_capacity (nel); \
01935           /* The optimal break-point as estimated from simulations */ \
01936           /* Note that Mergesort is O(nz log(nz)) while searching all */ \
01937           /* values is O(nr), where nz here is non-zero per row of */ \
01938           /* length nr. The test itself was then derived from the */ \
01939           /* simulation with random square matrices and the observation */ \
01940           /* of the number of non-zero elements in the output matrix */ \
01941           /* it was found that the breakpoints were */ \
01942           /*   nr: 500  1000  2000  5000 10000 */ \
01943           /*   nz:   6    25    97   585  2202 */ \
01944           /* The below is a simplication of the 'polyfit'-ed parameters */ \
01945           /* to these breakpoints */ \
01946           octave_idx_type n_per_col = (a_nc > 43000 ? 43000 : \
01947                                         (a_nc * a_nc) / 43000); \
01948           octave_idx_type ii = 0; \
01949           octave_idx_type *ri = retval.xridx(); \
01950           octave_sort<octave_idx_type> sort; \
01951           \
01952           for (octave_idx_type i = 0; i < a_nc ; i++) \
01953             { \
01954               if (retval.xcidx(i+1) - retval.xcidx(i) > n_per_col) \
01955                 { \
01956                   for (octave_idx_type j = a.cidx(i); j < a.cidx(i+1); j++) \
01957                     { \
01958                       octave_idx_type col = a.ridx(j); \
01959                       EL_TYPE tmpval = a.data(j); \
01960                       for (octave_idx_type k = m.cidx(col) ; \
01961                            k < m.cidx(col+1); k++) \
01962                         { \
01963                           octave_quit (); \
01964                           octave_idx_type row = m.ridx(k); \
01965                           if (w[row] < i + 1) \
01966                             { \
01967                               w[row] = i + 1; \
01968                               Xcol[row] = tmpval * m.data(k); \
01969                             } \
01970                           else \
01971                             Xcol[row] += tmpval * m.data(k); \
01972                         } \
01973                     } \
01974                   for (octave_idx_type k = 0; k < nr; k++) \
01975                     if (w[k] == i + 1) \
01976                       { \
01977                         retval.xdata(ii) = Xcol[k]; \
01978                         retval.xridx(ii++) = k; \
01979                       } \
01980                 } \
01981               else \
01982                 { \
01983                   for (octave_idx_type j = a.cidx(i); j < a.cidx(i+1); j++) \
01984                     { \
01985                       octave_idx_type col = a.ridx(j); \
01986                       EL_TYPE tmpval = a.data(j); \
01987                       for (octave_idx_type k = m.cidx(col) ; \
01988                           k < m.cidx(col+1); k++) \
01989                         { \
01990                           octave_quit (); \
01991                           octave_idx_type row = m.ridx(k); \
01992                           if (w[row] < i + 1) \
01993                             { \
01994                               w[row] = i + 1; \
01995                               retval.xridx(ii++) = row;\
01996                               Xcol[row] = tmpval * m.data(k); \
01997                             } \
01998                           else \
01999                             Xcol[row] += tmpval * m.data(k); \
02000                         } \
02001                     } \
02002                   sort.sort (ri + retval.xcidx(i), ii - retval.xcidx(i)); \
02003                   for (octave_idx_type k = retval.xcidx(i); k < ii; k++) \
02004                     retval.xdata(k) = Xcol[retval.xridx(k)]; \
02005                 }  \
02006             } \
02007           retval.maybe_compress (true);\
02008           return retval; \
02009         } \
02010     }
02011 
02012 #define SPARSE_FULL_MUL( RET_TYPE, EL_TYPE, ZERO ) \
02013   octave_idx_type nr = m.rows (); \
02014   octave_idx_type nc = m.cols (); \
02015   \
02016   octave_idx_type a_nr = a.rows (); \
02017   octave_idx_type a_nc = a.cols (); \
02018   \
02019   if (nr == 1 && nc == 1) \
02020     { \
02021       RET_TYPE retval = m.elem (0,0) * a; \
02022       return retval; \
02023     } \
02024   else if (nc != a_nr) \
02025     { \
02026       gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc); \
02027       return RET_TYPE (); \
02028     } \
02029   else \
02030     { \
02031       RET_TYPE retval (nr, a_nc, ZERO); \
02032       \
02033       for (octave_idx_type i = 0; i < a_nc ; i++) \
02034         { \
02035           for (octave_idx_type j = 0; j < a_nr; j++) \
02036             { \
02037               octave_quit (); \
02038               \
02039               EL_TYPE tmpval = a.elem(j,i); \
02040               for (octave_idx_type k = m.cidx(j) ; k < m.cidx(j+1); k++) \
02041                 retval.elem (m.ridx(k),i) += tmpval * m.data(k); \
02042             } \
02043         } \
02044       return retval; \
02045     }
02046 
02047 #define SPARSE_FULL_TRANS_MUL( RET_TYPE, EL_TYPE, ZERO, CONJ_OP ) \
02048   octave_idx_type nr = m.rows (); \
02049   octave_idx_type nc = m.cols (); \
02050   \
02051   octave_idx_type a_nr = a.rows (); \
02052   octave_idx_type a_nc = a.cols (); \
02053   \
02054   if (nr == 1 && nc == 1) \
02055     { \
02056       RET_TYPE retval = CONJ_OP (m.elem(0,0)) * a; \
02057       return retval; \
02058     } \
02059   else if (nr != a_nr) \
02060     { \
02061       gripe_nonconformant ("operator *", nc, nr, a_nr, a_nc); \
02062       return RET_TYPE (); \
02063     } \
02064   else \
02065     { \
02066       RET_TYPE retval (nc, a_nc); \
02067       \
02068       for (octave_idx_type i = 0; i < a_nc ; i++) \
02069         { \
02070           for (octave_idx_type j = 0; j < nc; j++) \
02071             { \
02072               octave_quit (); \
02073               \
02074               EL_TYPE acc = ZERO; \
02075               for (octave_idx_type k = m.cidx(j) ; k < m.cidx(j+1); k++) \
02076                 acc += a.elem (m.ridx(k),i) * CONJ_OP (m.data(k)); \
02077               retval.xelem (j,i) = acc; \
02078             } \
02079         } \
02080       return retval; \
02081     }
02082 
02083 #define FULL_SPARSE_MUL( RET_TYPE, EL_TYPE, ZERO ) \
02084   octave_idx_type nr = m.rows (); \
02085   octave_idx_type nc = m.cols (); \
02086   \
02087   octave_idx_type a_nr = a.rows (); \
02088   octave_idx_type a_nc = a.cols (); \
02089   \
02090   if (a_nr == 1 && a_nc == 1) \
02091     { \
02092       RET_TYPE retval = m * a.elem (0,0); \
02093       return retval; \
02094     } \
02095   else if (nc != a_nr) \
02096     { \
02097       gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc); \
02098       return RET_TYPE (); \
02099     } \
02100   else \
02101     { \
02102       RET_TYPE retval (nr, a_nc, ZERO); \
02103       \
02104       for (octave_idx_type i = 0; i < a_nc ; i++) \
02105         { \
02106           octave_quit (); \
02107           for (octave_idx_type j = a.cidx(i); j < a.cidx(i+1); j++) \
02108             { \
02109               octave_idx_type col = a.ridx(j); \
02110               EL_TYPE tmpval = a.data(j); \
02111               \
02112               for (octave_idx_type k = 0 ; k < nr; k++) \
02113                 retval.xelem (k,i) += tmpval * m.elem(k,col); \
02114             } \
02115         } \
02116       return retval; \
02117     }
02118 
02119 #define FULL_SPARSE_MUL_TRANS( RET_TYPE, EL_TYPE, ZERO, CONJ_OP ) \
02120   octave_idx_type nr = m.rows (); \
02121   octave_idx_type nc = m.cols (); \
02122   \
02123   octave_idx_type a_nr = a.rows (); \
02124   octave_idx_type a_nc = a.cols (); \
02125   \
02126   if (a_nr == 1 && a_nc == 1) \
02127     { \
02128       RET_TYPE retval = m * CONJ_OP (a.elem(0,0)); \
02129       return retval; \
02130     } \
02131   else if (nc != a_nc) \
02132     { \
02133       gripe_nonconformant ("operator *", nr, nc, a_nc, a_nr); \
02134       return RET_TYPE (); \
02135     } \
02136   else \
02137     { \
02138       RET_TYPE retval (nr, a_nr, ZERO); \
02139       \
02140       for (octave_idx_type i = 0; i < a_nc ; i++) \
02141         { \
02142           octave_quit (); \
02143           for (octave_idx_type j = a.cidx(i); j < a.cidx(i+1); j++) \
02144             { \
02145               octave_idx_type col = a.ridx(j); \
02146               EL_TYPE tmpval = CONJ_OP (a.data(j)); \
02147               for (octave_idx_type k = 0 ; k < nr; k++) \
02148                 retval.xelem (k,col) += tmpval * m.elem(k,i); \
02149             } \
02150         } \
02151       return retval; \
02152     }
02153 
02154 #endif
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Thu Mar 7 2013 16:03:21 for GNU Octave by  doxygen 1.7.1