Sparse-op-defs.h

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

Sparse-op-defs.h
