Mercurial > octave-nkf

comparison liboctave/Sparse-diag-op-defs.h @ 8964:f4f4d65faaa0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement sparse * diagonal and diagonal * sparse operations, double-prec only. Date: Sun, 8 Mar 2009 16:28:18 -0400 These preserve sparsity, so eye(5) * sprand (5, 5, .2) is *sparse* and not dense. This may affect people who use multiplication by eye() rather than full(). The liboctave routines do *not* check if arguments are scalars in disguise. There is a type problem with checking at that level. I suspect we want diag * "sparse scalar" to stay diagonal, but we have to return a sparse matrix at the liboctave. Rather than worrying about that in liboctave, we cope with it when binding to Octave and return the correct higher-level type. The implementation is in Sparse-diag-op-defs.h rather than Sparse-op-defs.h to limit recompilation. And the implementations are templates rather than macros to produce better compiler errors and debugging information.
author Jason Riedy <jason@acm.org>
date Mon, 09 Mar 2009 17:49:13 -0400
parents
children 1bba53c0a38d
comparison
equal deleted inserted replaced
8963:d1eab3ddb02d 8964:f4f4d65faaa0
1 /* -*- C++ -*-
2
3 Copyright (C) 2009 Jason Riedy, Jaroslav Hajek
4
5 This file is part of Octave.
6
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20
21 */
22
23 #if !defined (octave_sparse_diag_op_defs_h)
24 #define octave_sparse_diag_op_defs_h 1
25
26 template <typename RT, typename DM, typename SM>
27 RT do_mul_dm_sm (const DM& d, const SM& a)
28 {
29 const octave_idx_type nr = d.rows ();
30 const octave_idx_type nc = d.cols ();
31
32 const octave_idx_type a_nr = a.rows ();
33 const octave_idx_type a_nc = a.cols ();
34
35 if (nc != a_nr)
36 {
37 gripe_nonconformant ("operator *", nr, nc, a_nr, a_nc);
38 return RT ();
39 }
40 else
41 {
42 RT r (nr, a_nc, a.nnz ());
43
44 octave_idx_type l = 0;
45
46 for (octave_idx_type j = 0; j < a_nc; j++)
47 {
48 r.xcidx (j) = l;
49 const octave_idx_type colend = a.cidx (j+1);
50 for (octave_idx_type k = a.cidx (j); k < colend; k++)
51 {
52 const octave_idx_type i = a.ridx (k);
53 if (i >= nr) break;
54 r.xdata (l) = d.dgelem (i) * a.data (k);
55 r.xridx (l) = i;
56 l++;
57 }
58 }
59
60 r.xcidx (a_nc) = l;
61
62 r.maybe_compress (true);
63 return r;
64 }
65 }
66
67 template <typename RT, typename SM, typename DM>
68 RT do_mul_sm_dm (const SM& a, const DM& d)
69 {
70 const octave_idx_type nr = d.rows ();
71 const octave_idx_type nc = d.cols ();
72
73 const octave_idx_type a_nr = a.rows ();
74 const octave_idx_type a_nc = a.cols ();
75
76 if (nr != a_nc)
77 {
78 gripe_nonconformant ("operator *", a_nr, a_nc, nr, nc);
79 return RT ();
80 }
81 else
82 {
83
84 const octave_idx_type mnc = nc < a_nc ? nc: a_nc;
85 RT r (a_nr, nc, a.cidx (mnc));
86
87 for (octave_idx_type j = 0; j < mnc; ++j)
88 {
89 const typename DM::element_type s = d.dgelem (j);
90 const octave_idx_type colend = a.cidx (j+1);
91 r.xcidx (j) = a.cidx (j);
92 for (octave_idx_type k = a.cidx (j); k < colend; ++k)
93 {
94 r.xdata (k) = s * a.data (k);
95 r.xridx (k) = a.ridx (k);
96 }
97 }
98 for (octave_idx_type j = mnc; j <= nc; ++j)
99 r.xcidx (j) = a.cidx (mnc);
100
101 r.maybe_compress (true);
102 return r;
103 }
104 }
105
106 #endif // octave_sparse_diag_op_defs_h