GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
Sparse-perm-op-defs.h
Go to the documentation of this file.
1 /* -*- C++ -*-
2 
3 Copyright (C) 2009-2018 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
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License 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 <https://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 #include "octave-config.h"
27 
28 #include "PermMatrix.h"
29 #include "lo-array-errwarn.h"
30 #include "oct-locbuf.h"
31 #include "oct-sort.h"
32 #include "quit.h"
33 
34 // Matrix multiplication
35 
36 template <typename SM>
37 SM octinternal_do_mul_colpm_sm (const octave_idx_type *pcol, const SM& a)
38 // Relabel the rows according to pcol.
39 {
40  const octave_idx_type nr = a.rows ();
41  const octave_idx_type nc = a.cols ();
42  const octave_idx_type nent = a.nnz ();
43  SM r (nr, nc, nent);
44 
46 
47  for (octave_idx_type j = 0; j <= nc; ++j)
48  r.xcidx (j) = a.cidx (j);
49 
50  for (octave_idx_type j = 0; j < nc; j++)
51  {
52  octave_quit ();
53 
54  OCTAVE_LOCAL_BUFFER (octave_idx_type, sidx, r.xcidx (j+1) - r.xcidx (j));
55  for (octave_idx_type i = r.xcidx (j), ii = 0; i < r.xcidx (j+1); i++)
56  {
57  sidx[ii++]=i;
58  r.xridx (i) = pcol[a.ridx (i)];
59  }
60  sort.sort (r.xridx () + r.xcidx (j), sidx, r.xcidx (j+1) - r.xcidx (j));
61  for (octave_idx_type i = r.xcidx (j), ii = 0; i < r.xcidx (j+1); i++)
62  r.xdata (i) = a.data (sidx[ii++]);
63  }
64 
65  return r;
66 }
67 
68 template <typename SM>
69 SM octinternal_do_mul_pm_sm (const PermMatrix& p, const SM& a)
70 {
71  const octave_idx_type nr = a.rows ();
72  if (p.cols () != nr)
73  octave::err_nonconformant ("operator *",
74  p.rows (), p.cols (), a.rows (), a.cols ());
75 
76  return octinternal_do_mul_colpm_sm (p.col_perm_vec ().data (), a);
77 }
78 
79 template <typename SM>
80 SM octinternal_do_mul_sm_rowpm (const SM& a, const octave_idx_type *prow)
81 // For a row permutation, iterate across the source a and stuff the
82 // results into the correct destination column in r.
83 {
84  const octave_idx_type nr = a.rows ();
85  const octave_idx_type nc = a.cols ();
86  const octave_idx_type nent = a.nnz ();
87  SM r (nr, nc, nent);
88 
89  for (octave_idx_type j_src = 0; j_src < nc; ++j_src)
90  r.xcidx (prow[j_src]) = a.cidx (j_src+1) - a.cidx (j_src);
91  octave_idx_type k = 0;
92  for (octave_idx_type j = 0; j < nc; ++j)
93  {
94  const octave_idx_type tmp = r.xcidx (j);
95  r.xcidx (j) = k;
96  k += tmp;
97  }
98  r.xcidx (nc) = nent;
99 
100  octave_idx_type k_src = 0;
101  for (octave_idx_type j_src = 0; j_src < nc; ++j_src)
102  {
103  octave_quit ();
104  const octave_idx_type j = prow[j_src];
105  const octave_idx_type kend_src = a.cidx (j_src + 1);
106  for (k = r.xcidx (j); k_src < kend_src; ++k, ++k_src)
107  {
108  r.xridx (k) = a.ridx (k_src);
109  r.xdata (k) = a.data (k_src);
110  }
111  }
112  assert (k_src == nent);
113 
114  return r;
115 }
116 
117 template <typename SM>
118 SM octinternal_do_mul_sm_colpm (const SM& a, const octave_idx_type *pcol)
119 // For a column permutation, iterate across the destination r and pull
120 // data from the correct column of a.
121 {
122  const octave_idx_type nr = a.rows ();
123  const octave_idx_type nc = a.cols ();
124  const octave_idx_type nent = a.nnz ();
125  SM r (nr, nc, nent);
126 
127  for (octave_idx_type j = 0; j < nc; ++j)
128  {
129  const octave_idx_type j_src = pcol[j];
130  r.xcidx (j+1) = r.xcidx (j) + (a.cidx (j_src+1) - a.cidx (j_src));
131  }
132  assert (r.xcidx (nc) == nent);
133 
134  octave_idx_type k = 0;
135  for (octave_idx_type j = 0; j < nc; ++j)
136  {
137  octave_quit ();
138  const octave_idx_type j_src = pcol[j];
139  octave_idx_type k_src;
140  const octave_idx_type kend_src = a.cidx (j_src + 1);
141  for (k_src = a.cidx (j_src); k_src < kend_src; ++k_src, ++k)
142  {
143  r.xridx (k) = a.ridx (k_src);
144  r.xdata (k) = a.data (k_src);
145  }
146  }
147  assert (k == nent);
148 
149  return r;
150 }
151 
152 template <typename SM>
153 SM octinternal_do_mul_sm_pm (const SM& a, const PermMatrix& p)
154 {
155  const octave_idx_type nc = a.cols ();
156  if (p.rows () != nc)
157  octave::err_nonconformant ("operator *",
158  a.rows (), a.cols (), p.rows (), p.cols ());
159 
160  return octinternal_do_mul_sm_colpm (a, p.col_perm_vec ().data ());
161 }
162 
163 #endif
SM octinternal_do_mul_colpm_sm(const octave_idx_type *pcol, const SM &a)
SM octinternal_do_mul_pm_sm(const PermMatrix &p, const SM &a)
for large enough k
Definition: lu.cc:617
Sparse< bool > sort(octave_idx_type dim=0, sortmode mode=ASCENDING) const
Definition: Sparse.cc:2231
SM octinternal_do_mul_sm_pm(const SM &a, const PermMatrix &p)
SM octinternal_do_mul_sm_colpm(const SM &a, const octave_idx_type *pcol)
SM octinternal_do_mul_sm_rowpm(const SM &a, const octave_idx_type *prow)
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:400
double tmp
Definition: data.cc:6252
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
p
Definition: lu.cc:138
T * xdata(void)
Definition: Sparse.h:488
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:41
for i
Definition: data.cc:5264
octave_idx_type * xcidx(void)
Definition: Sparse.h:514