Mercurial > octave-nkf

comparison liboctave/Sparse-perm-op-defs.h @ 8968:91d53dc37f79

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add perm * sparse, perm \ sparse, sparse * perm, and sparse / perm operations. Nothing terribly fancy in any of this. There probably is some mechanism for using the permutation vectors and some assign or index method in the sparse classes, but I've never understood all the intricacies. I'm opting for a simple implementation at the cost of possibly duplicating some functionality.
author Jason Riedy <jason@acm.org>
date Tue, 10 Mar 2009 21:54:44 -0400
parents
children 829e69ec3110
comparison
equal deleted inserted replaced
8967:5bbbf482909a 8968:91d53dc37f79
1 /* -*- C++ -*-
2
3 Copyright (C) 2009 Jason Riedy
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_perm_op_defs_h)
24 #define octave_sparse_perm_op_defs_h 1
25
26 // Matrix multiplication
27
28 template <typename SM>
29 SM octinternal_do_mul_colpm_sm (const octave_idx_type *pcol, const SM& a)
30 // Relabel the rows according to pcol.
31 {
32 const octave_idx_type nr = a.rows ();
33 const octave_idx_type nc = a.cols ();
34 const octave_idx_type nent = a.nnz ();
35 SM r (nr, nc, nent);
36
37 for (octave_idx_type k = 0; k < nent; ++k)
38 {
39 OCTAVE_QUIT;
40 r.xridx (k) = pcol[a.ridx (k)];
41 r.xdata (k) = a.data (k);
42 }
43 for (octave_idx_type j = 0; j <= nc; ++j)
44 r.xcidx (j) = a.cidx (j);
45
46 r.maybe_compress (false);
47 return r;
48 }
49
50 template <typename SM>
51 SM octinternal_do_mul_pm_sm (const PermMatrix& p, const SM& a)
52 {
53 const octave_idx_type nr = a.rows ();
54 if (p.cols () != nr)
55 {
56 gripe_nonconformant ("operator *", p.rows (), p.cols (), a.rows (), a.cols ());
57 return SM ();
58 }
59
60 if (p.is_row_perm ())
61 {
62 // Form the column permutation and then call the colpm_sm routine.
63 const octave_idx_type *prow = p.pvec ().data ();
64 OCTAVE_LOCAL_BUFFER(octave_idx_type, pcol, nr);
65 for (octave_idx_type i = 0; i < nr; ++i)
66 pcol[prow[i]] = i;
67 return octinternal_do_mul_colpm_sm (pcol, a);
68 }
69 else
70 return octinternal_do_mul_colpm_sm (p.pvec ().data (), a);
71 }
72
73 template <typename SM>
74 SM octinternal_do_mul_sm_rowpm (const SM& a, const octave_idx_type *prow)
75 // For a row permutation, iterate across the source a and stuff the
76 // results into the correct destination column in r.
77 {
78 const octave_idx_type nr = a.rows ();
79 const octave_idx_type nc = a.cols ();
80 const octave_idx_type nent = a.nnz ();
81 SM r (nr, nc, nent);
82
83 for (octave_idx_type j_src = 0; j_src < nc; ++j_src)
84 r.xcidx (prow[j_src]) = a.cidx (j_src+1) - a.cidx (j_src);
85 octave_idx_type k = 0;
86 for (octave_idx_type j = 0; j < nc; ++j)
87 {
88 const octave_idx_type tmp = r.xcidx (j);
89 r.xcidx (j) = k;
90 k += tmp;
91 }
92 r.xcidx (nc) = nent;
93
94 octave_idx_type k_src = 0;
95 for (octave_idx_type j_src = 0; j_src < nc; ++j_src)
96 {
97 OCTAVE_QUIT;
98 const octave_idx_type j = prow[j_src];
99 const octave_idx_type kend_src = a.cidx (j_src + 1);
100 for (k = r.xcidx (j); k_src < kend_src; ++k, ++k_src)
101 {
102 r.xridx (k) = a.ridx (k_src);
103 r.xdata (k) = a.data (k_src);
104 }
105 }
106 assert (k_src == nent);
107
108 r.maybe_compress (false);
109 return r;
110 }
111
112 template <typename SM>
113 SM octinternal_do_mul_sm_colpm (const SM& a, const octave_idx_type *pcol)
114 // For a column permutation, iterate across the destination r and pull
115 // data from the correct column of a.
116 {
117 const octave_idx_type nr = a.rows ();
118 const octave_idx_type nc = a.cols ();
119 const octave_idx_type nent = a.nnz ();
120 SM r (nr, nc, nent);
121
122 for (octave_idx_type j = 0; j < nc; ++j)
123 {
124 const octave_idx_type j_src = pcol[j];
125 r.xcidx (j+1) = r.xcidx (j) + (a.cidx (j_src+1) - a.cidx (j_src));
126 }
127 assert (r.xcidx (nc) == nent);
128
129 octave_idx_type k = 0;
130 for (octave_idx_type j = 0; j < nc; ++j)
131 {
132 OCTAVE_QUIT;
133 const octave_idx_type j_src = pcol[j];
134 octave_idx_type k_src;
135 const octave_idx_type kend_src = a.cidx (j_src + 1);
136 for (k_src = a.cidx (j_src); k_src < kend_src; ++k_src, ++k)
137 {
138 r.xridx (k) = a.ridx (k_src);
139 r.xdata (k) = a.data (k_src);
140 }
141 }
142 assert (k == nent);
143
144 r.maybe_compress (false);
145 return r;
146 }
147
148 template <typename SM>
149 SM octinternal_do_mul_sm_pm (const SM& a, const PermMatrix& p)
150 {
151 const octave_idx_type nc = a.cols ();
152 if (p.rows () != nc)
153 {
154 gripe_nonconformant ("operator *", a.rows (), a.cols (), p.rows (), p.cols ());
155 return SM ();
156 }
157
158 if (p.is_row_perm ())
159 return octinternal_do_mul_sm_rowpm (a, p.pvec ().data ());
160 else
161 return octinternal_do_mul_sm_colpm (a, p.pvec ().data ());
162 }
163
164 #endif // octave_sparse_perm_op_defs_h